Christophe Louembet
email: louembet@laas.frprarante@laas.fr

Paulo R Arantes Gilz

Event-triggered Model predictive control for spacecraft rendezvous

Keywords: Event-triggered control, Predictive control, Impulsive control

In the context of space rendezvous hovering phase, a event-triggered predictive controller has been designed in order to maintain the relative position of a chaser spacecraft relatively to a target for orbit servicing purposes. The main goal is to minimize the consumption and the computation burden of the controller with respect to previously developed controller.

INTRODUCTION

Mastering the spacecraft rendezvous and its automatic control will open venue of economical opportunities for space industry. For instance, for orbit servicing operations, the ability of a space tug or supply spacecraft to maintain safely and thriftily its relative position with respect to the target spacecraft is crucial for different purposes (observations, waiting for order, etc). In this study, the target spacecraft is inert and the chaser one is moved by means of chemical thrusters so that the control is modelled by an impulsive signal. Such hovering capacities can be obtained by developing efficient control algorithms. Pursuing this aim, the authors have developed an impulsive predictive controller for steering to and maintaining the chaser spacecraft in the given polytopic subset Arantes [START_REF] Arantes Gilz | Model predictive control for rendezvous hovering phases based on a novel description of constrained trajectories[END_REF]. The control strategy accounts for the periodic nature of the relative motion between spacecraft to steer the chaser to the set of periodic orbits that are included in a particular polytope. In fact, in Arantes Gilz et al. (2018), this predictive impulsive controller has been proved to stabilize this set even in presence of saturations: this property has been obtained by considering at least three consecutive impulsive controls for sake of controllability. That MPC scheme induces a periodic computation and execution of the multi-impulse controls. However, if this predictive controller shows efficiency in numerous cases and ability to be embedded in space-certified computation board (LEON 3), it suffers few drawbacks inherent to its nature. For instance, since computed at a given period, the obtained controls can be unnecessary or too small to be executed by the thrusters especially if the system belongs or is closed to the admissible set. Moreover, since the convergence to the admissible set is only ensured after the three impulsive controls are applied, the invariance of the admissible set can not be set before the third impulse. To overcome those drawbacks, we propose to combine the periodically triggered controller presented above with the event-triggered controller exposed in this work. The role of the first controller is to steer the system close enough to the admissible set so that the proposed event-triggered controller can take over to stabilize the admissible set. The proposed controller will compute one impulse that P. Arantes Gilz is supported by the ANR FASTRelax.

will systematically bring the system to the admissible set when the suitable event occurs.

To develop this event-triggered predictive controller, several challenges have to be faced. First, The relevant signals to be observed must be defined and constructed. Trigger rules have also to be set carefully to ensure the good properties of invariance and avoid bad behavior such as Zeno phenomenon. Then, the computation of the impulse is addressed by means of optimization to ensure the convergence. Note that the trigger rules should lead to the feasibility of the mathematical program. Finally the cohesion between both controller has to be carefully studied in order to obtain a global convergence and local invariance of the admissible set. Reader must keep in mind that these goals has to be achieved while taking into account the periodic nature of the relative dynamics and the non trivial description of the admissible set. If this the abstract exposes our different choice of methods, several technical aspects of this work have been omitted for sake of brevity.

RENDEZVOUS CONTEXT, DYNAMIC MODEL AND TECHNOLOGICAL CONSTRAINTS

The relative motion consists in the motion of a chaser spacecraft equipped with thrusters with respect to the moving local frame attached to a passive target spacecraft. Under Keplerian assumptions, the relative motion between two spacecraft in the Earth gravitational field has been expressed by means of linearized impulsive differential equations:

Ẋ(t) = A(t)X(t) + B∆V i δ(t -t i ) with X(t) ∈ R 6
and A(t) being 2π-periodic. After a change of independent variable from time t to true anomaly ν of the target spacecraft e.g. its position on its orbit, a similar transformation is applied to obtain a second state space where the state are the parameters D exploited in several previous works [START_REF] Deaconu | Designing continuously constrained spacecraft relative trajectories for proximity operations[END_REF]; Arantes [START_REF] Arantes Gilz | Model predictive control for rendezvous hovering phases based on a novel description of constrained trajectories[END_REF]:

D (ν) = A D (ν)X(ν)+B D (ν)∆V i δ(ν -ν i )
. This state can be viewed as the coordinates of the relative orbits (see [START_REF] Deaconu | On the trajectory design, guidance and control for spacecraft rendezvous and proximity operations[END_REF], Chap 2) for full details). As stated earlier, the aim of the controller is to stabilize the set of relative periodic orbits satisfying the following relative position constraints:

x ≤ x(t) ≤ x, y ≤ y(t) ≤ y, z ≤ z(t) ≤ z , ∀t ≥ t 0 . This set is described by linear but time varying conditions on state D (see details in Arantes Gilz et al. ( 2017)):

S p D :=    D ∈ R 6 d 0 = 0, x ≤ M x (ν)D ≤ x y ≤ M y (ν)D ≤ y z ≤ M z (ν)D ≤ z , ∀ν    (1)
On top of that, the presence of a deadzone and saturations in the thrusters operation must be accounted for.

3. EVENT-BASED ALGORITHM Control Law For a given true anomaly ν, D + (ν), the state right after an impulse ∆V (ν) ∈ R 3 is given by D + (ν) = D(ν) + B D (ν)∆V (ν). We opt for the strategy where every control impulse produce systematically an periodic orbit such that:

d + 0 = d 0 + B D0 (ν)∆V = 0
(2) Equivalently, the control that satisfy equation ( 2) is given by ∆V = B ⊥ D0 (ν)λ + ∆V 0 (ν), λ ∈ R and B ⊥ D0 describes the kernel space of B D0 and ∆V 0 is a particular solution of (2). Assuming the periodicity pursuit strategy, the effect of a control impulse on the cuurent state D is described by

D + (ν, λ) = D(ν) + B D B ⊥ D0 (ν)λ + ∆V 0
(3) To maintain the state D in the admissible set S p D , the impulse control is computed by solving the following program: min

λ∈R (B ⊥ D0 (ν)λ) s.t. D + ∈ S p D λ ∈ I sat (Psat)
where I sat describes the input saturation and deadzone condition in function of λ such that

I sat (ν) = λ ∈ Rs.t. ∆V ≤ |B ⊥ D0 (ν)λ + ∆V 0 (ν)| ≤ ∆V (4)
Trigger laws The trigger law is designed to complete a threefold objective. First the problem (Psat) has to be feasible when called. Second, unnecessary controls must be avoided. Third Zeno phenomenon should be ensured to not occur. To set the trigger rules we need to define few terms. First, let the set ∆ + be set of state D reachable with one impulse and describe by (3). Since λ ∈ R, ∆ + is a timevarying line. A necessary condition for the admissible set to be reachable is that the line ∆ + intersects S p D . Let us denote this intersection Λ such that Λ = ∆ + ∩ S p D . S p D being convex bounded and closed (see Arantes Gilz et al. ( 2017)), Λ is a segment of the line ∆ + . If L(Λ), the length of segment Λ, is different from zero, then S p D is reachable without input constraints. Considering these constraints, the admissible set is reachable if the set Λ sat , defined by Λ sat = Λ ∩ I sat , is non empty. Note that Zeno behavior is avoided by the presence of a deadzone condition in (4). Using these definitions, trigger rules can be set: To resume, the application of the single impulse control obtained by solving (Psat) is triggered if the set Λ sat is vanishing i.e. when its length goes under a given thres-hold δ. Otherwise, if the admissible set is unreachable, the three impulse control can be computed as in Arantes [START_REF] Arantes Gilz | Model predictive control for rendezvous hovering phases based on a novel description of constrained trajectories[END_REF] and applied.

if D(ν) ∈ S p D : Wait. if D(ν) ∈ S

FIRST RESULTS AND DISCUSSION

First results are provided on figure 1. Scenario 1 of the work Arantes Gilz et al. ( 2018) is rerun with the combination of the multi-impulse controller and the proposed event-triggered single impulse one. The threshold δ is 0.005 and bounds ∆V and ∆V are resp. equal to 0.5cm • s -1 and 50cm • s -1 . For this scenario, the first impulses are computed by means of Arantes Gilz et al. ( 2018) controller. Then, the event triggered controller takes over at the end of the first period. From there, three other impulses are computed and applied each time the state lays outside of the admissible set while being reachable (L(Λ sat ) = 0). Then, when the system is close enough to the admissible set (L(Λ sat ) = 0), the single impulse controller takes over and stabilize the admissible set. The usage of such eventbased strategy permits to save 25% of fuel consumption with respect to the periodic multi-impulse approach. Several technical aspects have been omitted in this extended abstract. Mainly, the evaluation of the set Λ sat that is a part of the computation burden has not been exposed. In addition, if the stability of the combined controllers has not been formerly set, first hints at stability comes from the properties of the periodic multi-impulse controller.

  p D and Λ sat = ∅ Apply the three impulse strategy from Arantes Gilz et al. (2017) if D(ν) ∈ S p D , Λ sat = ∅, L(Λ sat ) < δ and d dν L(Λ sat ) < 0 Compute and apply ∆V = arg min((Psat)) else Wait.

Fig. 1 .

 1 Fig. 1. Observed signals, trigger and controlled trajectories


Cédric Maigrot

Vincent Claveau

Ewa Kijak

Ronan Sicre

Ronan

MediaEval 2016: A multimodal system for the Verifying Multimedia Use task

published or not. The documents may come

MediaEval 2016: A multimodal system for the Verifying Multimedia Use task

Conclusion

» Text-based approach: competes with the source-based approach in terms of recall but tends to classify every tweet as fake » Image-based approach: low precision compared with estimations on the training set.

This may be due to: (1) small and unbalanced reference database; (2) original image and forged ones are sometimes very similar; (3) presence of stamps » Combination-based approach: does not offer any gain due to overfitting

  Text-based approach (run-T) Detect if the message is style-wise similar to known hoax § Capture similar comments between an unknown image and an image from the training set (e.g. It's photoshopped ) and similar genres of comments (e.g. presence of smileys) § Prediction made by a k -Nearest-Neighbor approach (in this case k = 1) Source-based approach (run-S) Detect if the message is related to a trustworthy source § 2 type of sources searched: news-related organisms (e.g. press agencies) and explicit citations of the source of the image (e.g. the pattern photographed by + Name) § Predict real if a trustworthy source is detected, fake else Example Image-based approach (run-I) Detect a known image » Compare an unknown image to an image database of 8 000 known images (7 500 fake and 500 real images) » Database images extracted from 5 specialized websites » Description of an image by a deep CNN layer output (4096-dimensional descriptor) » Predict real (resp. fake) if a real (resp. fake) similar image is found in the database, uncertain else Combination approach (run-C) Combine the three previous predictions » Late fusion: learn the best combination » Boosting algorithm (adaboost.MH, parameters of the machine learning algorithm are set by cross-validation on the training data)

  Cédric Maigrot Vincent Claveau Ewa Kijak Ronan Sicre {firstname}.{lastname}@irisa.fr MediaEval 2016: A multimodal system for the Verifying Multimedia Use task

	Cédric Maigrot Vincent Claveau Ewa Kijak Ronan Sicre
		{firstname}.{lastname}@irisa.fr
	Why use a multimodal system ? Because there are several types of hoax !
	» False information present in the text content	» Forged image	» Image reused for an other event
	Global Hypotheses
	» Prediction is first made at the image-level, then propagated to the tweets that contain the image
	» Translation if the detected language is different than english

Acknowledgements

This work is partly supported by the Direction Générale de l'Armement, France (DGA).


Christophe Tirel

Marie-Charlotte C Renoult

Christophe Dumouchel

Viscoelastic

Viscoelastic Ligament Dynamics in Free Liquid Jet Experiments

Keywords: dilute polymer solution, ligament behavior, jet break-up

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

Introduction

Viscoelastic liquids have specific shear thinning and strain-hardening properties that confer on them characteristic behaviors while atomizing. For instance, laminar viscoelastic liquid jets break up with formation of beads-onstrings structures composed of large spherical beads and thin cylindrical threads. These threads might be the seat of an elongational flow where capillary and elastic forces are in balance, the first being the driving force, the second the resistive one, that leads to a specific decrease of their diameter with time. This decrease is controlled by the stress relaxation time scale tr of the viscoelastic solution. The knowledge of this time scale is therefore of paramount importance in the context of atomization but its experimental determination remains challenging especially for very dilute polymer solutions. The Capillary Breakup Extensional Rheometer (CaBER) performs this measurement by detecting and analyzing the elasto-capillary thinning of a filament formed by stretching a liquid bridge between two coaxial cylindrical plates. For dilute polymer solutions, Rodd and al. [1] have shown that the CaBER have limitations: there is a lower shear viscosity limit and a lower relaxation time limit for which the relaxation time measurement is not possible due to the fast thinning of the filament and the lack of resolution of the device. Moreover, liquid shape oscillations are introduced due to the fast transition of the plates from their initial positions to their final configurations for which the liquid bridge is unstable and thins down [2]. These oscillations persist for multiples of the capillary time and corrupt the CaBER measurement if the breakup occurs before they vanish. These limitations have been overcome in part with the use of the Slow Retraction Method and the use of a high speed camera, but there is a progressive appearance of beads-on-strings structures at the final steps of the filament thinning resulting from an iterated elastic instability [3]. Another technic based on jet experiments have been developed [4][5][6][START_REF] Bazilevsky | Third European Rheology Conference and Golden Jubilee Meeting of the British Society of Rheology[END_REF][START_REF] Keshavarz | [END_REF]. For dilute polymer solutions, the jet rheometer is seen as a promising alternative to the liquid filament rheometer [START_REF] Bazilevsky | Third European Rheology Conference and Golden Jubilee Meeting of the British Society of Rheology[END_REF]. Up so far, the jet rheometer considers the behavior of a cylindrical liquid jet forced at a given perturbation amplitude and frequency [4][5][6][START_REF] Keshavarz | [END_REF]. However, this technic remains tricky in particular because the resulting relaxation time might report a dependence on the jet operating conditions, i.e. the jet velocity and the initial disturbance amplitude and frequency [6]. Recently, new results suggest however the existence of an operating domain free of these dependences [START_REF] Keshavarz | [END_REF].

The present work intends to shed more light on the use of liquid jet capillary instability as a support to perform relaxation time measurement. Two characteristic features distinguish the present approach: 1the use of free liquid jets; 2the application of a multi-scale tool for analyzing the jets.

During the capillary instability of a viscoelastic cylindrical jet, the neck diameter Dmin(t) is expected to decrease according to three regimes that are presented in Table 1. The first regime is the Plateau-Rayleigh (PR) instability during which a sinusoidal perturbation develops and Dmin(t) decreases according to Eq. (1) where Djet is the jet diameter,  the non-dimensional growth rate and t the capillary time (= (Djet 3 /) 1/2  being the liquid density and  the surface tension) [9]. The second regime is the Elasto-Capillary (EC) regime during which elastic and capillary forces balance in the liquid strings whose diameter Dmin(t) exhibits an exponential decrease as reported by Eq. ( 2). The third regime is the Visco-Capillary (VC) regime during which elongational viscous and capillary forces balance in the liquid string. This mechanism is associated with a linear decrease of the neck diameter Dmin(t) (see Eq. ( 3) where tBU designates the jet breakup time and µe the terminal elongational viscosity. This presentation shows that the measurement of the relaxation time tr is possible provided that 1the EC regime exists; 2the temporal resolution is high enough to provide a exploitable Dmin(t) in this regime.

Table 1 Successive regimes of a viscoelastic jet capillary instability

Based on the behaviors summarized in Table 1, previous jet rheometers rested on the measurement of Dmin(t) which required the determination of its location. The present work uses a different approach through the application of a multi-scale analysis as introduced in previous papers [10,[START_REF] Tirel | Ilass Europe 28th european conference on liquid atomization and spray systems[END_REF]. This analysis describes the jet as a time-dependent scale distribution e2(d,t) and the analysis consists in considering the behavior of a small scale whose dynamics is the same as the diameter Dmin(t). The analysis of this dynamics allows identifying the existence of the EC regime and, then, measuring the solution relaxation time.

Compared to the previous works [10,[START_REF] Tirel | Ilass Europe 28th european conference on liquid atomization and spray systems[END_REF], the optical arrangement in the present work has been modified in order to improve the spatial resolution. It is indeed of paramount importance to have, at each time, enough scales (that is enough pixels) to correctly describe the ligaments that are the structures of interest in this study. Furthermore, by assuming the axisymmetry of the jets, the multi-scale analysis is performed in 3D which returns the temporal evolution of the specific area (amount of surface area per unit volume) of the jet. Finally, new injection conditions are considered here. More precisely, three discharge orifices with different diameters have been used. This article is structured as followed. The experimental setup and optical imaging arrangement are presented in the next section. It is followed by a section dedicated to the analyzing tools (image analysis and multi-scale analysis). Then, the results are presented and discussed. The article ends with a conclusion section.

Experimental setup

The experimental setup produces free falling cylindrical jets. The injector is composed of a cylindrical body (5 mm internal diameter) ended with a cylindrical discharge orifice of length Lor = 300 µm. Three different discharge orifice diameters are available, i. The choice of these resolutions results from a compromise between the number of pixels in the jets and the depth of field of the optical arrangement. The physical field covered by the image is 8.29 mm x 5.52 mm and 3.68 mm x 2.46 mm according to the spatial resolution. For each operating condition (orifice diameter, flow rate, position) n = 150 uncorrelated images are recorded. This number is high enough for a statistical analysis as detailed in the analysis section. Figure 1 presents typical images of the jets at the nozzle exit and in the breakup region for each orifice diameter.

In the near nozzle region, the white line in the middle of the jets attests that the flows are laminar and axisymmetric.

The appearance of the naturally-selected sinusoidal perturbation downstream from the nozzle is always clearly visible as well as the ligaments of the beads-on-strings pattern down to the breakup region.

Analysis

Viscoelastic flows issuing from a cylindrical orifice might be larger or smaller than the orifice diameter according to the characteristics of the internal flow. This is the reason why it is decided to characterize the jets with a proper measurement of their diameter Djet and of their subsequent velocity Vjet. For each injection flow rate, Djet is deduced from the jet projected surface area Sjet averaged on the 150 images and measured along a distance zjet = 3 mm where the jet is not deformed, i.e., Djet = Sjet/zjet. The corresponding velocity is obtained by: Vjet = 4Qv/Djet 2 . Furthermore, the wavelength jet of the naturally-selected mode identified in Fig. 1 2 Characteristics of the jets for all operating conditions (For case n° 3.1, the wavelength measurement was not possible)

As mentioned in the introduction, a 3D multi-scale tool is used to analyze the temporal evolution of the jets. This tool has been introduced in detail in previous articles [START_REF] Tirel | Ilass Europe 28th european conference on liquid atomization and spray systems[END_REF][START_REF] Dumouchel | [END_REF] and is summarized here only. It provides the scale distribution e3(d, z) which is the ratio of the surface area of the local system eroded at scale d on twice the total volume of this system before any erosion operation. This distribution is measured for all possible values of the scale d (i.e., until the jet has totally disappeared by erosion) and at several locations z from the nozzle exit to the breakup region. Although images provide a 2D information, the 3D scale distribution is deduced by assuming the axisymmetry of the liquid jet at every stages. The capillary instability is known to be axisymmetric enough for this assumption to be acceptable. (This is confirmed by the images in Fig. 1.) The steps of the scale distribution measurement are illustrated in Fig. 2. The initial image (Fig. 2-a) is treated in order to produce a two-level image (Fig. 2-b) where the liquid appears in black on a white background. Then, the Euclidian Distance Map of the jet is performed (Fig. 2-c). The volume-based scale distribution function is calculated in a portion of the jet delimited by an Analyzing Window (Fig. 2-d). d,z) represents the measurement performed on image i, the average scale distribution at this location is defined by:

𝑒 3 (𝑑, 𝑧) = ∑ 𝑉 𝑖 (𝑧) 𝑛 𝑖=1 𝑒 3𝑖 (𝑑, 𝑧) ∑ 𝑉 𝑖 (𝑧) 𝑛 𝑖=1 (4)
where Vi(z) is the jet portion volume in the Analyzing Window in image i. Then, the scale distribution spatial dependence is transformed as a temporal dependence by using the equivalent time t defined by:

𝑡 = √𝑉 𝑗𝑒𝑡 2 + 2𝑔𝑧 -𝑉 𝑗𝑒𝑡 𝑔 ( 5
)
where g is the gravitational acceleration. When Vjet 2 >> 2gz (Froude number >> 1), Eq. ( 5) simplifies as t = z/Vjet. The temporal resolution of the scale distribution measurement is depending on the shift of the Analyzing Window. For most operating conditions, using a shift equal to Hw appeared acceptable. However, when a greater temporal resolution is required (such as for injector Dor = 200 µm with small injection velocities) it was divided by 2, returning twice more points. The convergence of the mean scale distribution measurement, as a function of the image number, is studied through the value of e3(d=0,t), which is the maximum value of e3 by definition. The variation of the average of e3(0,t) is always less than 4% when n = 150 images (this result was established for condition 1.2, 1.3, 1.5, 1.6, 2.3, 2.5, 2.6, 2.7 and for at least three different positions for every case). n = 150 is therefore a sufficient number of images.

It is important to mention that the scale distribution measured as explained above describes the mean arrangement of the liquid system. This mean arrangement is called the unit pattern of the jet. If the jet is exactly cylindrical, the unit pattern is a portion of this cylinder whatever its length (e3(d) is independent of cylinder's length). As shown in previous works ( [START_REF] Dumouchel | ASME 2017 Fluids Engineering Division Summer Meeting[END_REF] for instance) e3(d) is linear in this case and decreases from 2/Dc at d = 0 to 0 at d = dmax = Dc where Dc stands for the cylinder diameter. If the cylindrical jet is deformed by a sinusoidal perturbation of wavelength , the unit pattern is a portion of the jet of length . Finally, if the jet is a succession of beads connected by cylindrical ligaments, the unit pattern is made of one bead whose diameter is equal to the mean diameter of the beads seen at this location and one ligament whose diameter is equal to the mean diameter of the ligaments seen at this location. Furthermore, the bead-ligament volume ratio of the unit pattern corresponds to the mean bead-ligament volume ratio seen at this location. An example of the temporal evolution of the scale distribution is given in Fig. 3 (Case n°1.6, refer to Table 2). At the initial time (t = 2.9 ms) e3(d,t) is linear indicating that the initial jet is cylindrical. This, of course, corresponds to the observation. As time goes, the maximum scale dmax (smallest scale for which e3(d,t) = 0) increases in agreement with the growth of the sinusoidal perturbation in Regime 1 (see Table 1) and the behavior of the beads in Regimes 2 and 3. For the great times, note that the scale distribution remains linear in the small region. This behavior denotes the persistence of a cylindrical portion in the unit pattern which, of course, corresponds to the string of the beads-on-strings structure. Two points are interesting to mention here: 1the derivative of the linear portion decreases with time which reveals the decrease of the ligament diameter (this is precisely this dynamic that we are interested in), 2 -e3(0,t) decreases with time. This second behavior is better shown in Fig. 4 for three cases (1.6; 2.6; 3.4). As mentioned earlier, e3(0,t) is half of the specific surface area (i.e. interface surface area per unit volume [START_REF] Evers | [END_REF]). As expected for a capillary instability this quantity decreases with time. This behavior is reported for every case investigated in this work. (The oscillations of e3(0,t) sometimes observed at the beginning, such as for case n°1.6 in Fig. 4, are likely related to the presence of a static deformation.)

Figure 3 Temporal evolution of e3(d,t) (Case n° 1.6) Figure 4 Temporal evolution of the specific surface e3(0,t) (Cases n° 1.6 ; 2.6 ; 3.4, the breakup times tBU are given in Table 2)

Now, the analysis consists in identifying scales whose dynamics is the same as the one of the unit pattern ligaments that are the seat of an elongational flow. Using the quasi one-dimensional continuity equation provided by Stelter et al. [15], the stretching rate 𝜀̇ of such a flow is equal to:

𝜀̇= 𝜕𝑣 𝜕𝑧 = -2 𝐷 𝑐 (𝑡) d𝐷 𝑐 (𝑡) d𝑡 ( 6
)
where Dc is the diameter of the cylindrical ligament and v the longitudinal velocity along the z axis. It is straight forward to demonstrate that for an isolated cylinder, the stretching rate given by Eq. ( 6) is related to the scale distribution by [START_REF] Dumouchel | ASME 2017 Fluids Engineering Division Summer Meeting[END_REF]:

𝜀̇= 𝑒 3 ′ ̇(𝑑, 𝑡) 𝑒 3 ′(𝑑, 𝑡) (7)
where the prime indicates a scale derivative and the dot a temporal derivative. Since the scale distribution is linear in the scale space for a cylinder, the stretching rate characterizing its elongation is independent of the scale d also.

For case 1.6, Fig. 5 shows the function 𝑒 3 ′ ̇(𝑑, 𝑡) at time t = 4.4 ms for which the jet is rearranged as beads-onstrings. This figure reveals that 𝑒 3 ′ ̇(𝑑, 𝑡) is scale-independent in the small scale region, i.e., the ligament thinning is identifiable with the scale distribution temporal evolution. Furthermore, the dynamics of this flow is the same as the dynamics of the scale d1 defined in Fig. 5 because of the scale proximity. Therefore, as done in previous works [10][START_REF] Tirel | Ilass Europe 28th european conference on liquid atomization and spray systems[END_REF][START_REF] Dumouchel | [END_REF][START_REF] Dumouchel | ASME 2017 Fluids Engineering Division Summer Meeting[END_REF], the temporal evolution of the necks of the jet can be explored from d1(t). However, in the present situation, an accurate determination of the scale d1 would require a spatial resolution that is difficult to achieve when imaging the full jet. The present optical arrangement provides d1(t) measurements with non-negligible scatter (see Fig. 6 for instance) that is pejorative for the analysis. This is the reason why it is decided to consider the scale d4 (see Fig. 5 for the definition) that is less sensitive to spatial resolution and whose dynamics remains the same as the scale d1 as evidenced in Fig. 6.

Results and Discussion

For six cases among the 21 that constitute the experimental part of this work (see Table 2), Fig 7 presents the temporal evolution of the characteristic scale d4(t) in two different manners: (d4(t)/Djet) and (1 -d4(t)/Djet). As expected, the ratio d4(t)/Djet highlights a monotonous decrease of the scale d4. For every operating condition, the curves in Fig. 7 are used to identify the Regime 1 (PR) and the Regime 3 (VC) introduced in Table 1. (Cases n° 1.2; 1.6; 2.2; 2.6; 3.2; 3.6)

According to Eq. ( 1), Regime 1 is associated with a linear increase of the function (1 -d4(t)/Djet) in a semi-log representation. Furthermore, being exponential (see Eq. ( 1)), this increase is expected to be the steepest one. This behavior is seen for every situation as evidenced by the dot lines in Fig. 7. To demonstrate that this linear increase corresponds to the PR regime, the theoretical increase expressed by Eq. ( 1) is calculated,  being the non-dimensional growth reported by Rayleigh linear theory [1] for the sinusoidal perturbation of wavelength . (The values of  are given in Table 2.) The continuous lines in Fig. 7 indicate the PR theoretical results. The experimental and theoretical growth rates are of the same order of magnitude which confirms that the identified linear portion does correspond to regime 1. (The results also exhibit a systematic smaller experimental growth rate compared to the theoretical one. It is interesting to mention that this very difference was reported by the analysis of simulated capillary instability of Newtonian ligaments [START_REF] Dumouchel | [END_REF].) The time t* at which regime 1 ends is identified as the time at which the linear behavior is lost (see Fig. 7).

According to Eq. ( 3), Regime 3 is associated with a linear decrease of the ratio (d4(t)/Djet) in a linear representation. Furthermore, this regime is the last one and should extend to the final breakup. The results presented in Fig. 7 show this behavior in every situation as indicated by the dash lines. The extension of this line allows identifying the time t** from which Regime 3 starts. (The time t** is shown in Fig. 7.)

The determination of the characteristic times t* and t** is performed for every operating condition and the results are displayed in Fig. 8. Following Table According to Eq. ( 2), the measurement of the relaxation time can be performed by considering the evolution of the scale d4(t) as a function of time in a semi-log representation. In order to compare the cases with each other, the ratio d4(t)/d4(t*) versus (tt*) is plotted in the interval [t*; t**] for the cases n° 1.3, 1.4, 1.5, 1.6. Figure 9 shows the results. Three important points must be mentioned. First, the four cases report a linear behavior in this representation which validates the dependence expressed by Eq. ( 2) and gives credit to the existence of the EC regime. Second, for every case, the number of points is sufficient to determine a slope with no ambiguity. Third, the slopes of the four cases are very close to each other. This means that the detected EC regime is not depending on the injection conditions and that the characteristic time that can be extracted from it is the relaxation time of the dilute polymer solution. For the present dilute solution, this relaxation time is tr = 0.71 ms.

Finally, Fig. 10 shows the ratio (t** -t*)/tr for every case. It can be seen that the cases reporting an EC regime are those for which the available temporal range [t*; t**] is at least 1.5 time the relaxation time tr. If it is less, this regime has no time to settle properly and to provide a possibility of measurement. On the other hand, we see that the jets whose breakup time are less than the relaxation time, no EC regime can be identified. This is the case of the smallest orifice diameter jets at every velocity for which the breakup time never exceeds 0.5 ms. The present experimental work and subsequent analysis demonstrate that a pure elongational flow is not guaranty in the strings of a beads-on-strings pattern formed during the breakup of free cylindrical-dilute-polymer-solution jets. One of the reason is that the successive regimes of small-scale dynamics exhibit short durations that sometimes are of the same order of magnitude as the characteristic time we want to measure or even less. Thus, the elongational flow has no time to properly settle and the corresponding elongational regime is unidentifiable. All this act as a brake upon the determination of dilute-polymer-solution relaxation time from the behavior of low-velocity cylindrical jets. These results might explain some dependences between the relaxation time tr and the injection conditions reported in the literature. Indeed, the regimes of the small-scale dynamics are recognized by their expected temporal variation. However, their small duration associated with small variations of the scale of interest might generate confusion between the tracked mathematical behaviors. In a previous work [START_REF] Tirel | Ilass Europe 28th european conference on liquid atomization and spray systems[END_REF], we suggested a possible link between a bead coalescence mechanism along the jet and the dependence between tr and the operating conditions. The simulation work due to Li and Fontelos [16] demonstrates that the coalescence mechanism doesn't modify the liquid strings' dynamics. Therefore, this mechanism is not responsible for the hypothetic presence of an EC regime. This, of course, doesn't exclude the fact that the absence of such a regime modifies the coalescence mechanism. A confident measurement of the dilute-polymer-solution relaxation time from a jet experiment is possible. This is what this work demonstrates. It is not clear yet which association of diameter and velocity should be recommended to perform this measurement, but the experimental work presented here might define the bases of a protocol to perform it with no ambiguity. The two important aspects of this protocol are 1the use of free liquid jets (meaning that it can be applied on any ligamentary structures whatever their production process); 2the use of a multiscale description tool that straight-forwardly applies in any conditions.

Regime 1 : 2 : 3 :

 123 Plateau-Rayleigh (PR) Regime Elasto-Capillary (EC) Regime Visco-Capillary (VC) 𝐷 𝑗𝑒𝑡 -𝐷 𝑚𝑖𝑛 (𝑡) ∝ 𝑒𝑥𝑝 (

Figure 1

 1 Figure 1 Viscoelastic jets at the nozzle exit and in the breakup region. a) Dor = 200 µm, Qv = 5.5 ml/min; b) Dor = 105 µm, Qv = 2.72 ml/min; c) Dor = 42 µm, Qv = 1.54 ml/min Transmission images of the jets are performed. In order to get a high spatial resolution, a 29 Mpixels camera (6576 x 4384 pixels, SVS-Vistek HR29050) is used. The spatial resolution has been settled to 1.26 µm/pixel for the

Figure 2

 2 Figure 2 Steps of scale distribution measurement a) Raw image; b) Two-level image; c) Euclidian Distance Map; d) Analyzing Window where the local scale distribution is measuredThe spatial evolution of the volume-based scale distribution is obtained by sliding the Analyzing Window from the nozzle exit down to the breakup region. The height Hw of the analyzing window is taken of the order of the naturally-selected perturbation wavelength. This choice is a good compromise between localness and number of images required to have a converged average scale distribution. Indeed, the local measurement of the scale distribution is performed on each image of the series and a local average distribution is calculated. Ife3i(d,z) represents the measurement performed on image i, the average scale distribution at this location is defined by:

Figure 5 Figure 6

 56 Figure 5 Temporal evolution of 𝑒 3 ′ ̇(d,t) (Case 1.6).Definition of characteristic scales

Figure 7

 7 Figure 7 Temporal evolution of scale d4. Identification of Regime 1 (PR) and Regime 3 (VC)(Cases n° 1.2; 1.6; 2.2; 2.6; 3.2; 3.6)

Figure 8

 8 shows that the range of this interval is very much dependent on the working conditions. For the small orifice diameter (Dor = 42 µm), the range of [t*; t**] is nil whatever the velocity indicating the absence of an EC regime for these jets. For the medium orifice diameter (Dor = 105 µm), the range of [t*; t**] remains pretty small (especially at high velocity) and might not allow the development of the EC regime. For the large orifice diameter (Dor = 200 µm), the range of [t*; t**] appears small for the two smallest and the highest velocities. However, for the four medium velocities (cases n° 1.3, 1.4, 1.5, 1.6), the range of [t*; t**] might be large enough to allow the development of the EC regime and, subsequently, the measurement of the relaxation time tr.

Figure 8 Figure 9

 89 Figure 8 Characteristic times t* and t** for every working condition

Figure 10

 10 Figure 10 Comparison of the [t*; t**] range with the solution relaxation time tr ConclusionThe present experimental work and subsequent analysis demonstrate that a pure elongational flow is not guaranty in the strings of a beads-on-strings pattern formed during the breakup of free cylindrical-dilute-polymer-solution jets. One of the reason is that the successive regimes of small-scale dynamics exhibit short durations that sometimes are of the same order of magnitude as the characteristic time we want to measure or even less. Thus, the elongational flow has no time to properly settle and the corresponding elongational regime is unidentifiable. All this act as a brake upon the determination of dilute-polymer-solution relaxation time from the behavior of low-velocity cylindrical jets. These results might explain some dependences between the relaxation time tr and the injection conditions reported in the literature. Indeed, the regimes of the small-scale dynamics are recognized by their expected temporal variation. However, their small duration associated with small variations of the scale of interest might generate confusion between the tracked mathematical behaviors. In a previous work[START_REF] Tirel | Ilass Europe 28th european conference on liquid atomization and spray systems[END_REF], we suggested a possible link between a bead coalescence mechanism along the jet and the dependence between tr and the operating conditions. The simulation work due to Li and Fontelos[16] demonstrates that the coalescence mechanism doesn't modify the liquid strings' dynamics. Therefore, this mechanism is not responsible for the hypothetic presence of an EC regime. This, of course, doesn't exclude the fact that the absence of such a regime modifies the coalescence mechanism. A confident measurement of the dilute-polymer-solution relaxation time from a jet experiment is possible. This is what this work demonstrates. It is not clear yet which association of diameter and velocity should be recommended to perform this measurement, but the experimental work presented here might define the bases of a protocol to perform it with no ambiguity. The two important aspects of this protocol are 1the use of free liquid jets (meaning that it can be applied on any ligamentary structures whatever their production process); 2the use of a multiscale description tool that straight-forwardly applies in any conditions.

  is measured and averaged on 50 images for each operating condition. All these measurements are reported in Table2.

	n°	Qv (ml/min)	Djet (µm)	Vjet (m/s)	tBU (ms)	tσ (ms)	λjet (µm)	ω (-)
				Dor = 200 µm
	1.1	3.56	195	1.99	11.1 0.137 978 0.336
	1.2	5.50	190	3.23	9.98 0.132 951 0.336
	1.3	7.96	185	4.94	7.62 0.127 975 0.330
	1.4	10.3	184	6.46	6.45 0.126 842 0.342
	1.5	12.2	200	6.49	6.88 0.142 975 0.339
	1.6	14.2	187	8.61	6.56 0.129 938 0.336
	1.7	17.3	190	10.2	5.02 0.132 911 0.340
				Dor = 105 µm
	2.1	1.75	125	2.37	6.03 0.070 578 0.342
	2.2	2.72	121	3.94	3.91 0.067 587 0.340
	2.3	3.93	117	6.09	3.21 0.064 569 0.337
	2.4	5.08	114	8.29	2.53 0.061 522 0.335
	2.5	7.02	112	11.9	2.06 0.060 506 0.333
	2.6	8.72	112	14.8	1.84 0.060 484 0.337
	2.7	10.2	112	17.3	1.73 0.060 480 0.340
				Dor = 42 µm
	3.1	1.09	46.0 10.9 0.571 0.016	-	-
	3.2	1.54	47.0 14.8 0.534 0.016 199 0.339
	3.3	1.75	47.0 16.8 0.518 0.016 198 0.338
	3.4	2.12	46.0 21.3 0.486 0.016 195 0.339
	3.5	2.36	47.0 22.7 0.457 0.016 187 0.328
	3.6	2.54	46.5 24.9 0.523 0.016 199 0.340
	3.7	2.72	46.5 26.7 0.506 0.016 202 0.341
	Table

  1, Regime 2i.e. the one of interest in the present workmay extend in the interval [t*; t**] at best.


Thomas Dewez
email: t.dewez@brgm.fr

Cécile Allanic

Gabriel Courrioux

Imed Ksibi

Simon Lopez

Gildas Noury

Silvain Yart

From underground laser scans to 3D urban geological and geotechnical models

Keywords: handheld laser scanner, underground quarries, collapse hazard, geomodelling, urban geological model, BIM

The near sub-surface geology, say down to 20-30-m-depth, of many cities has been massively exploited for extracting building stones and various other industrial or agricultural materials (gypsum, lime, etc…). The long-term instability of these cavities poses a significant collapse hazard conditioned by their geometry (void location, dimensions and shape) and by their surrounding rock mechanics properties. In this presentation, we show how handheld laser scanning surveys efficiently document geometric variables and can interact with 3D geological modelling of the surrounding rocks. The construction of near-surface urban geological models can then be turned into 3D geotechnical models by attributing geotechnical parameters to rock horizons and ultimately become a key subsurface knowledge component of BIM (Building Information Model).

Acquiring surface and subsurface geometry is no longer a challenge thanks to handheld laser scanners. Survey loop traverses can be pieced together to link surface and subsurface geometry with accuracies better than 1 m (an accuracy level compatible with urban risk management maps at 1/5.000) [START_REF] Dewez | Towards cavity collapse hazard maps with Zeb-Revo handheld laser scanner point clouds[END_REF]. However, the hundreds of-millions of 3D points describing the cavity surface cannot be integrated as such into geomodeling software. Too many points with not high enough information. We suggest two different scenarios to perform their integration: (i) as independent validation of geomodeling hypotheses, or (ii) as geomodel constraints.

In the first integration scenario, point cloud information is passed to the geomodeling software at a minimal level. A decimated triangular meshed model can be used to intersect the geomodel. Triangulation is performed at the point cloud processing software level (e.g. GeoSLAM desktop or Cloud Compare) and intersection is handled at the geomodelling software level with a generic query concept (here GeoModeller software with a generic query API -LOISELET et al.,2016). In this instance, cavity mesh triangular faces are refined based on the geological model queries (relying on the marching triangles algorithm) and provide geotechnical attributes based on the geological formations given by the geomodel. This scenario offers a visual display of geological properties (Fig. 1) for checking that modelled layers and structures match those observed in underground outcrops.

In the second scenario, which is more integrated, higher level information is passed to the geomodelling tool. Planar surfaces of marker horizons are segmented from the point cloud either manually using Compass (THIELE et al., 2017) or semi-automatically with FACETS [START_REF] Dewez | FACETS : a Cloud Compare plugin to extract geological planes from unstructured 3D point clouds[END_REF] and passed as structural data objects to constrain the geomodel (Fig. 2). This data integration is demonstrated on a ca. 1 ha underground building stone quarry of the eastern suburbs of Orléans, Central France. The cavity was scanned at ca. 1pt/1cm with a Zeb-Revo (90 Mpts underground and 35 Mpts above ground). A geomodel of the subsurface area (Calcaire de Beauce, Tertiary) was created with the GeoModeller software as a tabular sub-horizontal multilayer environment. The geomodel infers rock distribution over a domain of ca. 200 x 200 m with geological and geotechnical information (e.g. limit pressure for dimensioning building foundations).

Both approaches leverage a generic API query tool informing which domain surrounds a point and whether a geological contact cross-cuts a triangular face.

Figure 1 :

 1 Figure 1: Intersecting the geological model (GeoModeller) with the cavity wall 3D mesh. Modelled layers interpolate limit pressure values (p l *) spatially. Weathered : p l * 0.5 -1 MPa; Upper facies: p l * 1.5-2.5 MPa; Calcaire de Beauce : p l * 4.0 MPa


Martin Campos
email: campos@ljll.math.upmc.fr

Martin Campos Pinto

Moment preserving local spline projection operators

Keywords: Splines, Stable dual bases, Projection operators, Preservation of moments Mathematics Subject Classification (2010) 65D07, 41A15, 65D15

Moment preserving local spline projection operators

Introduction

Since the early developments of splines spaces and their description as linear combinations of locally supported B-splines {ϕ i : 0 ≤ i < n}, many local dual bases {ψ i : 0 ≤ i < n} have been proposed, see e.g. [START_REF] De Boor | The method of projections as applied to the numerical solution of two point boundary value problems using cubic splines[END_REF][START_REF] De Boor | On uniform approximation by splines[END_REF][START_REF] De Boor | Spline approximation by quasi-interpolants[END_REF][START_REF] De Boor | On local linear functionals which vanish at all B-splines but one[END_REF][START_REF] Lyche | Local spline approximation methods[END_REF][START_REF] Thomas | Bézier projection: A unified approach for local projection and quadrature-free refinement and coarsening of NURBS and T-splines with particular application to isogeometric design and analysis[END_REF][START_REF] Buffa | On Quasi-Interpolation Operators in Spline Spaces[END_REF]. Thanks to wellestablished stability properties these dual bases provide local projection operators P : f → i ψ i , f ϕ i that are stable in L 1 or L ∞ , and are also high order accurate thanks to the polynomial reproduction properties of the splines. In particular, they may be used as efficient alternatives to the non-local L 2 projection or nodal interpolation operators.

In some cases, it is important that the dual basis itself has high order polynomial reproduction properties: on an abstract level this guarantees that the dual projection Q : f → i ϕ i , f ψ i also satisfies some high order approximation properties, and for the primal spline projection P it corresponds to a moment preserving property. In mortar methods for instance [START_REF] Bernardi | Asymptotic and numerical methods for partial differential equations with critical parameters[END_REF], the use of local dual bases spanning high order polynomials allows to localize the matching conditions at the interfaces without loosing the optimality of the method, see [START_REF] Wohlmuth | A Mortar Finite Element Method Using Dual Spaces for the Lagrange Multiplier[END_REF][START_REF] Oswald | On polynomial reproduction of dual FE bases[END_REF][START_REF] Dornisch | Dual and approximate dual basis functions for B-splines and NURBS -Comparison and application for an efficient coupling of patches with the isogeometric mortar method[END_REF]. Another framework where such a tool is useful is the extension of classical commuting diagrams [START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF][START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF][START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF][START_REF] Buffa | Isogeometric discrete differential forms in three dimensions[END_REF] to non-conforming discretization spaces, see [START_REF] Pinto | Gauss-compatible Galerkin schemes for timedependent Maxwell equations[END_REF][START_REF] Pinto | Constructing exact sequences on non-conforming discrete spaces[END_REF][START_REF] Pinto | Compatible Maxwell solvers with particles I: conforming and non-conforming 2D schemes with a strong Ampere law[END_REF] and the presentation below.

One specific feature of B-splines compared to standard Finite Element spaces is their high order smoothness. This reduces their localization and makes it difficult to apply the abstract construction tools proposed in [START_REF] Oswald | On polynomial reproduction of dual FE bases[END_REF]. On the other hand, local dual bases with polynomial reproduction properties have been studied in the framework of spline wavelet multiresolutions [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF][START_REF] Dahmen | Biorthogonal Spline Wavelets on the Interval-Stability and Moment Conditions[END_REF] but their structure is made complicated by the refinability properties required by the multiresolution analysis, and their extension to non-uniform spline spaces is far from straightforward.

In this article we thus propose a direct and elementary construction of local dual bases for B-splines that reproduce polynomials of any desired degree. These dual bases are discontinuous but stable in L ∞ , and they share the same (possibly non-uniform) piecewise polynomial structure as the splines. As a consequence, the resulting local projection operators are moment preserving, stable and near optimal in L q , 1 ≤ q ≤ ∞, and they can be implemented with simple quadrature formulas for the spline coefficients.

The outline is as follows. In the remainder of Section 1 we review a few classical tools that we will need, and specify some notation. Section 2 then specifies the role played by moment preserving projection operators in the design of non-conforming commuting diagrams, which was an important motivation for this work. A couple of classical and not-so-classical local dual bases are then presented in Section 3, and our new dual basis is presented in Section 4, together with its stability analysis. Section 5 then provides a couple of numerical tests to verify the practical validity of the proposed construction.

Univariate B-splines

We begin by collecting some well-known facts about splines which can be found in [START_REF] De Boor | A practical guide to splines[END_REF]Ch. IX] or [START_REF] Schumaker | Spline functions: basic theory[END_REF]Ch. 4], and set some notation. Given two integers p, n, and a knot vector

ξ = {0 = ξ 0 ≤ • • • ≤ ξ n+p = 1}, (1)
we remind that B-splines B p i = B p [ξ i , . . . , ξ i+p+1 ] can be recursively defined on the interval Ω = [0, 1] by the Cox-de Boor formula, namely B 0 i = 1 [ξ i ,ξ i+1 ) for 0 ≤ i < n + p, and for q = 1, . . . p,

B q i (x) = x -ξ i ξ i+q -ξ i B q-1 i (x) + ξ i+q+1 -x ξ i+q+1 -ξ i+1 B q-1 i+1 (x), 0 ≤ i < n + p -q.
In this article we restrict ourselves to the case where the spline space

S p,ξ := Span{B p i } n-1 i=0 (2)
contains every p-degree polynomial and is of maximal regularity, which corresponds to the following assumption.

Assumption 1

The extremal knots have multiplicities p + 1 (the knot vector is said open), and the inner knots have multiplicities 1, i.e.

0 = ξ 0 = • • • = ξ p < ξ p+1 < • • • < ξ n-1 < ξ n = • • • = ξ n+p = 1.
Under Assumption 1 we may denote by

I k := [ζ k , ζ k+1 ) with ζ k := ξ k+p , 0 ≤ k < N := n -p, (3)
the non-empty intervals corresponding to the knot vector [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF]. The spline space then corresponds to

S p,ξ = {v ∈ C p-1 (Ω) : v| I k ∈ P p (I k ), 0 ≤ k < N }
where P p is the space of polynomials of degree less or equal to p, see [START_REF] Schumaker | Spline functions: basic theory[END_REF]Ex. 4.3], and the B-splines

ϕ i := B p i , 0 ≤ i < n, (4)
form a basis of S p,ξ with local supports supp(ϕ i ) = [ξ i , ξ i+p+1 ], see [START_REF] Schumaker | Spline functions: basic theory[END_REF]Th. 4.9].

Local projection operators with moment preservation properties

A convenient way of building local projection operators on S p,ξ that preserve polynomial moments up to some given degree m is to design functions

ψ i ∈ L ∞ (Ω), 0 ≤ i < n, (5)
with the following properties:

they are dual to the B-splines (4),

ψ i , ϕ j = δ i,j , 0 ≤ i, j < n, (P1)
they are local, in the sense that there exists K = K(p, m) such that supp ψ i ⊂ I + i with

I + i = I + i (K) := |k-i|≤K I k (P2)
they are stable, i.e. there exists a constant C = C(p, m) such that

ψ i L ∞ ≤ Ch -1 i where h i := |ξ i+p+1 -ξ i |, (P3)
they span every polynomial of degree m,

P m (Ω) ⊂ Span{ψ i } n-1 i=0 . (P4)
Thus, our problem is to design a dual basis with the properties listed above. To propose a solution we shall make the additional assumption that the grid is locally quasi-uniform.

Assumption 2 There exists a constant c * (possibly depending on p and m) such that

1 c * ≤ I k I k-1 ≤ c * , k = 1, . . . , N -1. ( 6
)
Using the above properties it is indeed standard to prove the following result, following for instance the arguments of Th. 6.60 in [START_REF] Schumaker | Spline functions: basic theory[END_REF].

Theorem 3 Let p, n, m ∈ N and let ξ be a knot vector satisfying Assumption 1 and 2. If (P1)-(P3) hold for given functions ψ i , 0 ≤ i < n, then the operator

P : L 1 (Ω) → S p,ξ , f → n-1 i=0 ψ i , f ϕ i (7)
is a projection on S p,ξ which is locally near-optimal in the sense that for 1 ≤ q ≤ ∞,

f -P f L q (I k ) ≤ C inf v∈S p,ξ f -v L q (I * k ) , 0 ≤ k < N (8)
holds with

I * k = ∪ k≤i≤k+p I + i ⊂ I + k (K + p)
, and a constant C that only depends on p and m. Moreover if (P4) holds, then P preserves the polynomial moments of degree m, g, P f = g, f ,

∀g ∈ P m , f ∈ L 1 (Ω). ( 9
)
Proof The projection property is straightforward to verify using (P1), and the moment preserving property simply follows by noticing that any g ∈ P m writes g = 0≤i<n γ i ψ i according to (P4), therefore

g, P f = 0≤i<n g, ϕ i ψ i , f = 0≤i<n γ i ψ i , f = g, f , f ∈ L 1 (Ω).
For [START_REF] De Boor | A practical guide to splines[END_REF] we first establish a local stability estimate. Given k, we observe that ϕ i vanishes on I k for i / ∈ {k, . . . , k + p}, and 0 ≤ ϕ i ≤ 1 holds by definition. Hence we have ϕ i L q (I k ) ≤ |I k | 1/q ≤ h 1/q i , and

P f L q (I k ) ≤ k+p i=k | ψ i , f | ϕ i L q (I k ) ≤ k+p i=k | ψ i , f |h 1 q i .
The products are then bound using a Hölder inequality, (P2) and (P3): it yields

| ψ i , f | ≤ ψ i L q (I + i ) f L q (I + i ) ≤ ψ i L ∞ |I + i | 1 q f L q (I + i ) ≤ Ch -1+ 1 q i f L q (I + i )
where we have used the local quasi-uniformity property ( 6) and (P2) to estimate |I + i | ≤ Ch i with a constant independent of i and N . As -1 + 1 q = -1 q , the above bounds show that

P f L q (I k ) ≤ C f L q (I * k ) (10)
holds with a constant C depending only on p. The local error estimate is then easily obtained by observing that

f -P f L q (I k ) ≤ f -v L q (I k ) + P (f -v) L q (I k ) ≤ (1 + C) f -v L q (I * k )
holds for all v ∈ S p,ξ .

Non-conforming commuting diagrams involving weak derivatives

An important motivation for moment preserving projection operators comes from the observation that they play a central role in commuting diagrams for nonconforming discretizations, which themselves are a key tool in the stable discretization of many systems of PDEs, see e.g. [START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF][START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF][START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF][START_REF] Buffa | Isogeometric discrete differential forms in three dimensions[END_REF][START_REF] Pinto | Gauss-compatible Galerkin schemes for timedependent Maxwell equations[END_REF].

To describe this setting with univariate spline discretizations, we denote by

V h = S p,ξ and W h = S p-1,ξ
two spline spaces of degree p and p-1 with open knot vectors ξ and ξ corresponding to the same grid ζ 0 , . . . , ζ N , under Assumption 1. This allows us to consider the standard derivative as a discrete operator mapping V h to W h , denoted

d h : V h → W h , v → ∂ x v
and its discrete adjoint d * h : W h → V h defined by the relations

d * h w, v = w, d h v v ∈ V h , w ∈ W h .
Formally the operator -d * h can be seen as a weak approximation of ∂ x on H 1 0 (Ω), and it is easy to verify that if one denotes by P V h the L 2 projection on V h , namely

P V h f, v = f, v , f ∈ L 2 (Ω), v ∈ V h , (11)
and similarly for P W h , then the following diagram commutes

H 1 0 (Ω) L 2 (Ω) W h V h ∂ x P W h P V h -d * h ( 12
)
in the sense that -d * h P W h u = P V h ∂ x u holds for all u ∈ H 1 0 (Ω). As described in [START_REF] Pinto | Gauss-compatible Galerkin schemes for timedependent Maxwell equations[END_REF][START_REF] Pinto | Constructing exact sequences on non-conforming discrete spaces[END_REF], commuting diagrams can be extended to non-conforming discretizations in order to gain some locality and flexibility in the design of the Finite Element spaces. Here a natural choice for the non-conforming space consists of broken splines. Given a division of the grid into subdomains containing each Ñ cells (assuming N = Ñ S for some S ∈ N) of the form

Ĩs = [ ζs , ζs+1 ] with ζs := ζ s Ñ , 0 ≤ s < S,
and writing ξs = { ξs i : 0 ≤ i ≤ ñ + p} with ñ = Ñ + p the open knot vector corresponding to p degree splines on the nodes

ζ k , s Ñ ≤ k ≤ (s + 1) Ñ , we let Ṽh := {v ∈ L 2 (Ω) : v| Ĩs ∈ S p, ξs , 0 ≤ s < S} (13)
denote the resulting broken spline space, which satisfies V h = Ṽh ∩ C p-1 (Ω). A projection P of the form ( 7) can then be used as a local smoothing operator mapping Ṽh on V h , and this allows to define a non-conforming differential operator dh : Ṽh → W h , ṽ → ∂ x P ṽ and its discrete adjoint d * h : W h → Ṽh characterized by d * h w, ṽ = w, ∂ x P ṽ , ṽ ∈ Ṽh .

A natural question is whether replacing d * h , V h and P V h by their non-conforming counterparts d * h , Ṽh and P Ṽh (the L 2 projection on Ṽh ) in ( 12) results in a commuting diagram. The answer is negative, and instead one must consider the operator

P * : L 2 (Ω) → Ṽh , P * f, ṽ = f, P ṽ , ṽ ∈ Ṽh . ( 14
)
Indeed for all u ∈ H 1 0 (Ω) and ṽ ∈ Ṽh , we have d * h P W h u, ṽ = -P W h u, ∂ x P ṽ = -u, ∂ x P ṽ = ∂ x u, P ṽ = P * ∂ x u, ṽ which means that the following diagram commutes:

H 1 0 (Ω) L 2 (Ω) W h Ṽh ∂ x P W h P * -d * h (15)
We note that applying P * amounts to applying P on the broken B-splines and to inverting mass matrices in the subdomains. As the latter have a bounded number of cells this is a local computation, just as applying the L 2 projection on Ṽh ,

P Ṽh f, ṽ = f, ṽ , ṽ ∈ Ṽh , (16)
amounts to solving local problems on the subdomains. We also observe that if Q : f → n-1 i=0 ϕ i , f ψ i denotes the dual projection operator mentioned in the introduction, we have P * = P Ṽh Q. Interestingly, P * cannot be a projection if Ṽh involves more than one subdomain: indeed that would imply ũ, ṽ = P * ũ, ṽ = ũ, P ṽ , ũ, ṽ ∈ Ṽh and taking ũ = ṽ = (I -P ) w for an arbitrary w ∈ Ṽh would give (I -P ) w = 0, hence w = P w ∈ V h , in contradiction with the fact that V h is a proper subspace of Ṽh . However, one easily verifies that if P is a moment preserving projection then P * is a quasi-interpolation.

Theorem 4 Let p, n, m ∈ N with m ≤ p, and let ξ be a knot vector satisfying Assumptions 1 and 2. Given S, Ñ ∈ N such that S Ñ = N := n -p, let then Ṽh be the broken spline space defined by [START_REF] Pinto | Gauss-compatible Galerkin schemes for timedependent Maxwell equations[END_REF]. If (P1)-(P4) hold for given functions ψ i , 0 ≤ i < n, then the operator P * : L 2 → Ṽh defined by ( 14) is a quasi-interpolation,

P * g = g, g ∈ P m , (17)
and it satisfies a local error estimate

f -P * f L q ( Ĩs ) ≤ C inf g∈P m f -g L q ( Ĩ * s ) , 0 ≤ s < S,
with a constant C that depends on p and Ñ the number of elements in a subdomain, and Ĩ * s ⊂ I + s Ñ ( K), K = K + max( Ñ , p), and interval which consists of a bounded number of cells close to the node ζ s Ñ , see (P2).

Proof To verify the quasi-interpolation property we use ( 14) and ( 9) to compute P * g, ṽ = g, P ṽ = g, ṽ , ṽ ∈ Ṽh .

The equality P * g = g then follows from the fact that for m ≤ p, P m is a subspace of Ṽh . Turning to the a priori estimate, we denote by { φs,a ∈ L 2 ( Ĩs ) : 0 ≤ a < ñ} an orthonormal basis for the local spline space S p, ξs , 0 ≤ s < S, and set φs,a = h 1 2 s 1 Ĩs φs,a and ψs,a = h-1 s φs,a with hs = | Ĩs |.

These functions satisfy ψs,a , φt,b = δ (s,a),(t,b) for 0 ≤ s, t < S, 0 ≤ a, b < ñ, and standard arguments based on the piecewise polynomial structure of the space S p, ξs over its quasi-uniform grid (6) allow us to write φs,a L q ≤ C h 1 q s and ψs,a

L q ≤ C h-1+ 1 q s ( 18
)
with constants depending only on p and Ñ , the number of polynomial pieces in each subdomain. By definition of P * we then have

P * f = S-1 s=0 ñ-1 a=0
P ψs,a , f φs,a .

A local stability estimate can then be obtained for P * using similar arguments as above: On a subdomain Ĩs we estimate

P * f L q ( Ĩs ) ≤ ñ-1 a=0 | P ψs,a , f | h 1 q s
and we infer from (P2) that the support of any P ψs,a = 0≤i<n ψ i , ψs,a ϕ i with 0 ≤ a < ñ consists of a bounded number of intervals around ζs = ζ s Ñ , that is Ĩ * s := supp(P ψs,a )

⊂ i :I + (K) i ∩ Ĩs =∅   i-p≤k<i I k   ⊂ I + s Ñ ( K)
with K = K + max( Ñ , p) as announced. It follows that we can bound

| P ψs,a , f | ≤ P ψs,a L q f L q ( Ĩ * s ) ≤ C h-1+ 1 q s f L q ( Ĩ * s )
where we have used the L q stability of P , namely Eq. ( 10), and the estimates [START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF]. Using again -1 + 1 q = -1 q , the above bounds give

P * f L q ( Ĩs ) ≤ C f L q ( Ĩ * s )
with a constant C depending only on p, q and Ñ . The local error estimate is then easily derived from the quasi-projection property [START_REF] Dornisch | Dual and approximate dual basis functions for B-splines and NURBS -Comparison and application for an efficient coupling of patches with the isogeometric mortar method[END_REF], indeed

f -P * f L q ( Ĩs ) ≤ f -g L q ( Ĩs ) + P * (f -g) L q ( Ĩs ) ≤ (1 + C) f -g L q ( Ĩ * s )
holds for all g ∈ P m .

Local dual bases

Before describing a solution to our problem we recall a couple of local dual bases, i.e. functions ψ i , 0 ≤ i < n, that satisfy properties (P1)-(P3) but not (P4).

The dual basis of de Boor

A classical L ∞ -stable local dual basis is described in [START_REF] De Boor | On local linear functionals which vanish at all B-splines but one[END_REF], see also [START_REF] Schumaker | Spline functions: basic theory[END_REF]Th. 4 , i = 0, 1, . . . , p + 1 (called perfect because its p-th derivative is of constant absolute value on its support [-1, 1]). A scaled transition function is then defined as

G i (x) = g 2x -ξ i -ξ i+p+1 ξ i+p+1 -ξ i with g(x) =      0 if x < -1, x -1 B p * if -1 ≤ x < 1, 1 if 1 ≤ x
and the resulting dual basis functions read

ψ db i (x) = d p+1 dx p+1 G i (x) i+p j=i+1 (x -ξ j ) p! , 0 ≤ i < n. (19)
It is known (see, e.g. the proof of Theorem 4.41 in [START_REF] Schumaker | Spline functions: basic theory[END_REF]) that these functions satisfy properties (P1) and (P3), and one easily verifies that supp(ψ db i ) ⊂ [ξ i , ξ i+p+1 ], so that (P2) holds with K = p.

Dual functions based on the piecewise polynomial structure of the B-splines

Another dual basis is easily derived from the piecewise polynomial structure of the B-splines. The resulting projection operator corresponds to the Bézier projection introduced in [START_REF] Thomas | Bézier projection: A unified approach for local projection and quadrature-free refinement and coarsening of NURBS and T-splines with particular application to isogeometric design and analysis[END_REF] and is studied as a particular case of local quasi-interpolation operators in [START_REF] Buffa | On Quasi-Interpolation Operators in Spline Spaces[END_REF].

An attractive feature of this basis is that it has the same polynomial structure as the B-splines, which simplifies the use of accurate quadrature formulas in numerical implementations. It can be defined as follows. For 0 ≤ k < N , denote

ϕ i,k := ϕ i | I k , 0 ≤ i < n, (20)
the restriction of the B-spline

ϕ i = B p i to the cell I k = [ζ k , ζ k+1 ).
For k ≤ i ≤ k +p these functions form a basis of P p (I k ), see e.g. [START_REF] Schumaker | Spline functions: basic theory[END_REF]Th. 4.21]: we may then let

ψ i,k ∈ P p (I k ), k ≤ i ≤ k + p, (21)
denote the dual polynomial basis on I k , characterized by the relations

ψ j,k , ϕ i,k I k = δ i,j for k ≤ i, j ≤ k + p.
The resulting dual functions are then defined as

ψ ps i := k∈K(i) α i,k 1 I k ψ i,k with α i,k := ϕ i -1 I k ϕ i (22)
where K(i) = k ∈ {0, . . . , N -1} : k ≤ i ≤ k + p , 0 ≤ i < n, denotes the cell indexes where the B-spline ϕ i has a non zero contribution.

Lemma 5 For 0 ≤ k < N , the mapping Π I k : f → k+p j=k ψ j,k , f I k ϕ j,k coincides with the L 2 (orthogonal) projection on P p (I k ).

Proof Writing ψ j,k = k+p l=k c k j,l ϕ l,k we see that

δ i,j = ψ j,k , ϕ i,k I k = k+p l=k c k j,l ϕ l,k , ϕ i,k I k ,
which also shows that the matrix (c k j,l ) k≤j,l≤k+p is symmetric. It follows that

k+p j=k ϕ i,k , ϕ j,k I k ψ j,k , f I k = k+p j,l=k ϕ i,k , ϕ j,k I k c k l,j ϕ l,k , f I k = ϕ i,k , f I k
which characterizes Π I k as the orthogonal projection on P p (I k ).

Remark 6

In particular, the linear forms λ I k , = 1, . . . , p + 1 from [10, Eq. ( 11)] and their analogs from [22, Eq. ( 58)] correspond to

λ I k (f ) = ψ i,k , f I k with i = k + -1.
Theorem 7 Let p ∈ N. The functions ψ ps i defined by [START_REF] Thomas | Bézier projection: A unified approach for local projection and quadrature-free refinement and coarsening of NURBS and T-splines with particular application to isogeometric design and analysis[END_REF] satisfy the properties (P1)-(P3), i.e. they form a local and stable dual basis to the B-splines. Moreover they satisfy (P4) with m = 0, i.e. they span the constant functions.

Proof By definition of the set K(i), we compute

k∈K(i) α i,k = k∈K(i) ϕ i -1 I k ϕ i = 1, (23)
hence

ψ ps i , ϕ j = k∈K(i) α i,k ψ i,k , ϕ j,k I k = k∈K(i) α i,k δ i,j = δ i,j
which proves the duality property (P1). The locality property (P2) is verified by construction. For the stability (P3) we can use Theorem 1 in [START_REF] Buffa | On Quasi-Interpolation Operators in Spline Spaces[END_REF] thanks to Remark 6 above. Under the local quasi-uniformity Assumption (2), this result states that

| ψ i,k , f I k | ≤ C|I k | -1 f L 1 (I k )
holds for k ≤ i ≤ k + p with a constant depending only on p. This allows to write

ψ i,k L ∞ = sup f ∈L 1 ψ i,k , f I k f L 1 (I k ) ≤ C|I k | -1 ≤ Ch -1 i ,
where the last inequality uses again the local quasi-uniformity [START_REF] De Boor | On local linear functionals which vanish at all B-splines but one[END_REF]. This clearly proves (P3) as the coefficients α i,k are all bounded (by 1). To finally prove the conservation property we observe that ϕ i α i,k = 1, ϕ i,k is the i-th coefficient of the function 1 in the local basis [START_REF] Schumaker | Spline functions: basic theory[END_REF]. In particular, we have

n-1 i=0 ϕ i ψ ps i = n-1 i=0 k∈K(i) 1, ϕ i,k 1 I k ψ i,k = N -1 k=0
1 I k = 1 Ω
which proves that (P4) indeed holds with m = 0.

Remark 8 As studied in [START_REF] Thomas | Bézier projection: A unified approach for local projection and quadrature-free refinement and coarsening of NURBS and T-splines with particular application to isogeometric design and analysis[END_REF] and [START_REF] Buffa | On Quasi-Interpolation Operators in Spline Spaces[END_REF], any bounded combination of the form [START_REF] Thomas | Bézier projection: A unified approach for local projection and quadrature-free refinement and coarsening of NURBS and T-splines with particular application to isogeometric design and analysis[END_REF] with k∈K(i) α i,k = 1 provides a stable and local dual basis for the B-splines. In [START_REF] Thomas | Bézier projection: A unified approach for local projection and quadrature-free refinement and coarsening of NURBS and T-splines with particular application to isogeometric design and analysis[END_REF] the values α i,k = ( ϕ i ) -1 I k ϕ i are advocated as leading to more accurate projections, and they are the only ones for which the dual basis spans the constants.

Dual functions based on a macro-subdivision

We now present our solution to the problem described in Section 1.2.

Construction of the dual basis

To preserve an arbitrary number of moments when projecting a function in the spline space, we construct a dual basis that uses a subdivision of the domain Ω into macro-elements containing each m + p + 1 cells of the form We next subdivide the B-splines ϕ i , 0 ≤ i < n, in two groups.

I k = [ζ k , ζ k+1 ), see (3)
-B-splines ϕ i that do not vanish on a macro-vertex ζ = ζ M will be said to be of macro-vertex type. Under Assumption 1 they are of the form

ϕ i with i ∈ I mv := { M, • • • , M + p -1} for some 0 ≤ < L.
-B-splines ϕ i that have their support contained in a single macro-element Î will be said to be of macro-element type. They are of the form

ϕ i with i ∈ I me := { M + p, • • • , ( + 1)M -1} for some 0 ≤ < L.
We gather the previous index sets into

I mv := 0≤ ≤L I mv and I me := 0≤ <L I me
and define the macro-subdivision dual basis functions ψ ms i as follows.

-If i ∈ I mv is the index of a macro-vertex spline, then we may take any stable dual basis function for the associated macro-vertex dual function. In order to have the same polynomial structure as the B-spline, we choose to set

ψ ms i := ψ ps i for i ∈ I mv . (25)
-If i ∈ I me is the index of a macro-element spline, we proceed in two steps. Given 0 ≤ < L, we let {φ j : j ∈ I me } be the polynomial basis of P m ( Î ) that satisfies φ j , ϕ l = δ j,l for j, l ∈ I me .

Here the important point is that such a basis is indeed well defined, and stable: it is guaranteed by Lemma 9 below. The macro-element dual function is then defined as

ψ ms i := φ i - j∈I mv ∪I mv +1 φ i , ϕ j ψ ms j for i ∈ I me , 0 ≤ < L. (27)
Lemma 9 Let 0 ≤ < L. There exists a unique basis {φ i : i ∈ I me } of P m ( Î ) that satisfies the relations

φ i , ϕ j = δ i,j i, j ∈ I me , (28)
moreover the estimate

φ i L ∞ ( Î ) ≤ Ch -1 i , i ∈ I me (29)
holds with a constant that only depends on p and m.

We postpone the proof of this lemma to Section 4.2 and state the desired result.

Theorem 10 Let p, m ∈ N. The functions ψ ms i defined by (25) and (27) satisfy the properties (P1)-(P4), i.e. they form a local and stable dual basis to the B-splines that span the polynomials of degree ≤ m.

Proof The duality property (P1) is straightforward to verify, given that the functions ψ ps i also satisfy this property. To verify the locality we first use the fact that the functions ψ ps i have the same local supports (ξ i , ξ i+p+1 ) as the B-splines, hence it is also the case for the dual functions of macro-vertex type. For the dual functions of macro-element type we observe that the functions φ i vanish by construction outside the macro-elements Î . It follows that

supp(ψ ms i ) ⊂ Î ∪ j∈I mv ∪I mv +1 supp(ψ ps j ) = (ξ M , ξ ( +1)M +2p )
which proves (P2) with K = 2p + m. The stability property (P3) is then easily derived from Theorem 7 and Lemma 9, using also the local quasi-uniformity (6) for the bound

| φ i , ϕ j | ≤ φ i L ∞ ϕ j L 1 ≤ Ch -1 i h j ≤ C
which holds with constants C, C depending only on p and m.

Stability analysis

We now prove Lemma 9. Given 0 ≤ < L we denote

φa = ϕ p+ M +a , 0 ≤ a ≤ m, (30)
the macro-element B-splines associated with Î = Î = [ ζ0 , ζM ]. We then let M ∈ M m+1 (R) be the matrix of moments

M a,b := 1 a!(m -a)! ζM ζ0 ( ζM -x) a ( ζM -ζ0 ) a φb (x) dx, 0 ≤ a, b ≤ m.
Verifying (28) amounts to showing that M is invertible: the polynomials φ i are then given by

φ p+ M +b = 0≤a≤m η b,a a!(m -a)! ( ζM -x) a ( ζM -ζ0 ) a
with η b,a := (M -1 ) b,a .

We will actually show that the inversion is stable, namely that

M -1 ≤ C| Î| -1 (31)
with a constant depending only on p and m. This in turn will allow us to write

φ p+ M +b L ∞ ≤ 0≤a≤m |η b,a |C m ≤ C| Î| -1
which leads to the desired bound (29), thanks to the local quasi-uniformity [START_REF] De Boor | On local linear functionals which vanish at all B-splines but one[END_REF].

To prove (31) we consider an arbitrary vector γ = (γ a ) 0≤a≤m ∈ R m+1 and let g := 0≤a≤m γ a φa be the associated spline. The vector µ := M γ then corresponds to the moments of g: its coefficients read

µ a = 1 a!(m -a)! ζM ζ0 ( ζM -x) a ( ζM -ζ0 ) a g(x) dx, 0 ≤ a ≤ m. (32)
We next let G 0 := g and define

G a+1 (x) := x ζ0
G a (y) dy for a = 0, . . . , m.

By construction, G := G m+1 is a spline of degree M on Î = [ ζ0 , ζM ], therefore it belongs to the space S M, ξ associated with the open knot vector ξ = { ξ0 , . . . , ξ3M } obtained by repeating M + 1 times the extremal nodes ζ0 and ζM , namely

ξ0 = • • • = ξM < • • • < ξ2M = • • • = ξ3M and ξM+k = ζk for 0 ≤ k ≤ M. (33)
It will be convenient to denote by Bq i the B-splines of degree q ≤ M on these knots. The macro-element B-splines (30) then correspond to

φa = Bp M +a , 0 ≤ a ≤ m, (34)
and since G is a linear combination of M -degree B-splines on these knots, we may write

G = 2M -1 i=0 Γ i BM i .
One can then verify that the first M coefficients vanish in this sum: indeed by construction one has

G (r) = G m+1-r for 0 ≤ r ≤ m g (r-m-1) for m + 1 ≤ r ≤ m + p = M -1, (35)
so that G (r) ( ζ0 ) = 0 for 0 ≤ b ≤ M -1. This guarantees that G only involves B-splines BM i with indexes i ≥ M , which can also be seen as a consequence of the de Boor-Fix formula [START_REF] De Boor | Spline approximation by quasi-interpolants[END_REF][START_REF] De Boor | Splines as linear combinations of B-splines. A survey[END_REF]

Γ i = M r=0 (-1) M -r Ψ (M -r) i (τ i )G (r) (τ i ) with Ψ i (x) = ( ξi+1 -x) • • • ( ξi+M -x) M ! .
(36) Here τ i can be arbitrary in ( ξi , ξi+M+1 ), and we further notice that any knot ξi+k with 1 ≤ k ≤ M is allowed (even if not in the open interval): indeed in the sum (36) the only discontinuous term corresponds to r = M , and on these knots the function Ψ (M -r) i = Ψ i vanishes, which makes the sum actually continuous. Thus, for 0 ≤ i < M one may take τ i = ξM = ζ0 and the formula yields Γ i = 0 indeed. For M ≤ i ≤ 2M -1 we then take τ i = ξ2M = ζM , which yields

|Γ i | ≤ M -1 r=0 | Î| r r! |G (r) ( ζM )|.
Using again that g consists of macro-element B-splines, the relations (35) imply that G (r) ( ζM ) = 0 if m + 1 ≤ r. For 0 ≤ r ≤ m we then compute

G (r) ( ζM ) = G m+1-r ( ζM ) = ζM ζ0 G m-r (x) dx = ζM ζ0 ( ζM -x)G m-r-1 (x) dx (if r ≤ m -1) . . . = ζM ζ0 ( ζM -x) m-r (m -r)! g(x) dx = r! | Î| m-r µ m-r
where the last equality follows from (32). This yields

Γ ∞ ≤ C m | Î| m µ . (37)
Turning to the function g, we use the recurrence relations (5.2) in [START_REF] De Boor | Splines as linear combinations of B-splines. A survey[END_REF] and the equality (34) to write

m a=0 γ a Bp M +a = g = G (m+1) = 2M -1 i=M Γ i BM i (m+1) = 2M -1 i=M Γ (m+1) i Bp i with differentiated coefficients defined by Γ (r) i := Γ i if r = 0 M +1-r ξi+M+1-r -ξi Γ (r-1) i -Γ (r-1) i-1 if r ≥ 1. (38)
This shows that γ a = Γ

M +a , and in particular Γ

M +a = 0 for a > m. Moreover, arguing as above one verifies that G (r) only involves splines BM-r i with i ≥ M , hence the recurrence relations [7, Eq. (5.2)] yield Γ (r) i = 0 for all i < M and 0 ≤ r ≤ m+1. According to (38) this allows us to only consider vectors of the form (Γ (r) M +a ) m a=0 , (which we shall conveniently denote by Γ (r) ) when we evaluate the differentiated coefficients Γ (m+1) M +a . This restriction is helpful since for 0 ≤ a ≤ m and 1 ≤ r ≤ m + 1, the interval ( ξM+a , ξ2M+a+1-r ) always contains at least one non-empty cell I k of the original grid. The local quasi-uniformity property (6) then allows us to write

|Γ (r) a | ≤ c| Î| -1 Γ (r-1)
for 1 ≤ r ≤ m + 1, with a constant c = c(p, m). Using (37) and (38) this leads to

γ = Γ (m+1) ≤ c m | Î| -m-1 Γ ≤ C| Î| -1 µ
with another constant C = C(p, m). This proves (31) and ends the proof of Lemma 9.

Numerical results

We end our study with a numerical assessment of the convergence properties of the various approximation operators described above.

Accuracy of the local projection operators P on the spline space

In Figure 1 we show the L 2 convergence curves of the local projection operators P on the smooth spline space V h = S p,ξ , as defined by [START_REF] De Boor | Splines as linear combinations of B-splines. A survey[END_REF]. The different operators tested here are reminded in Table 1 and their accuracy is compared to that of the L 2 projection P V h which is optimal. The plots on the left correspond to the approximation of the H 1 function f abs (x) = |x -1/π|, and the plots on the right correspond to the approximation of the smooth (C ∞ ) Runge function f runge (x) = 1/(25(2x -1) 2 + 1). Both are computed with degrees p ∈ {1, 2, 3}, from top to bottom.

Overall the accuracies of the different projection operators are very similar, and close to that of the optimal L 2 projection, which confirms the result of Theorem 3. For f abs (on the left) some differences are visible when cubic splines are employed: in particular the moment preserving P ms projection has a reduced accuracy on the coarsest meshes, which is due to its oscillating behavior close to the point where the target function is not C 1 . For the smooth f runge (on the right) the curves are hardly distinguishable Fig. 1 Convergence curves corresponding to the projection operators P : L 2 → V h defined on smooth spline spaces of specified degrees 1 ≤ p ≤ 3, as indicated. Definitions of the different operators are reminded in Table 1.

Accuracy of the commuting approximations P * on the broken splines

In Figure 2 we show the L 2 convergence curves of the operators P * which provide a commuting diagram for the broken spline space Ṽh , as studied in Section 2. The definitions of the different operators are reminded in Table 2 and their accuracy is compared to that of the L 2 projection P Ṽh which is optimal. Again we note that the latter operator is local because of the broken nature of the space, but it does not yield a commuting diagram, as discussed in Section 2. The target functions and spline degrees are the same as in Figure 1.

In contrast to Figure 1, we find that the accuracies of the different operators are not comparable, even for the smooth target function f runge . As announced by our a priori estimate from Theorem 4 the operator P * ms involving macro-subdivision Eq. ( 7), ( 19) ps local projection Pps based on the spline polynomial structure Eq. ( 7), ( 22) ms local projection Pms based on the macro-spline subdivision Eq. ( 7), ( 25), (27

)
dual functions ψ ms i seems to converge with an optimal rate, in accordance with their high order moment preserving properties. The operator P * ps involving ψ ps i dual functions has a lesser accuracy but still converges, in accordance with its first order moment preserving properties (see Theorem 7), and the operator P * db involving ψ db i dual functions hardly converges for degrees p ≥ 2.

Table 2 Approximation operators Pkey : L 2 → Ṽh used in Figure 2, for given values of key.

key operator reference L2 local L 2 projection P Ṽh Eq. ( 16) db local operator P * db using de Boor dual functions Eq. ( 14), ( 19) ps local operator P * ps using ps dual functions Eq. ( 14), ( 22) ms local operator P * ms using ms dual functions Eq. ( 14), ( 25), ( 27) Fig. 2 Convergence curves corresponding to the projection operators P * : L 2 → Ṽh defined on broken spline spaces with specified degrees.

  . Assuming that such a subdivision is possible, i.e., that there exists a positive integer L such that N = LM with M := m + p + 1 (under Assumption 1 we have N = n -p), we denote the elements of the macrosubdivision by Î = [ ζ , ζ +1 ) with ζ = ζ M . (24)

Table 1

 1 Projection operators P key : L 2 → V h used in Figure 1, for given values of key. key operator reference L2 non-local L 2 projection P V h Eq. (11) db local projection of de Boor P db

Acknowledgements The author would like to thank Eric Sonnendrücker and Ahmed Ratnani for the fruitful discussions that motivated this research.