Castro Carlos
email: carlos.castro@upm.es.

Nicolae Cindea
email: nicolae.cindea@math.univ-bpclermont.fr

Carlos Castro

Nicolae Cîndea

Arnaud Münch
email: arnaud.munch@math.univ-bpclermont.fr.

Controllability of the linear 1D wave equation with inner moving forces

Keywords: Linear wave equation, null controllability, finite element methods, Mixed formulation Mathematics Subject Classification (2010)-35L10, 65M12, 93B40

published or not. The documents may come

Controllability of the linear 1D-wave equation with inner moving forces

Introduction

Let T be a positive real, Q T the domain (0, 1) × (0, T ), q T a non-empty subset of Q T and Σ T = {0, 1} × (0, T ). We are concerned in this work with the null distributed controllability for the 1D wave equation:

   y tt -y xx = v 1 q T ,
(x, t) ∈ Q T y(x, t) = 0, (x, t) ∈ Σ T y(x, 0) = y 0 (x), y t (x, 0) = y 1 (x),

x ∈ (0, 1).

(

) 1
We assume that (y 0 , y 1 ) ∈ V := H 1 0 (0, 1) × L 2 (0, 1); v = v(t) is the control (a function in L 2 (q T )) and y = y(x, t) is the associated state. 1 q T from Q T to {0, 1} denotes the indicatrice function of q T . We also use the notation:

L y := y tty xx .

(

) 2
For any (y 0 , y 1 ) ∈ V and any v ∈ L 2 (q T ), there exists exactly one solution y to [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF], with the regularity y ∈ C 0 ([0, T ]; H 1 0 (0, 1)) ∩ C 1 ([0, T ]; L 2 (0, 1)) (see [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF]).

INTRODUCTION

The null controllability problem for (1) at time T is the following: for each (y 0 , y 1 ) ∈ V , find v ∈ L 2 (q T ) such that the corresponding solution to (1) satisfies y(• , T ) = 0, y t (• , T ) = 0 in (0, 1).

(3)

When the subset q T takes the form q T := ω × (0, T ), where ω denotes a subset of (0, 1), the null-controllability of (1) at any large time T > T ⋆ is well-known (for instance, see [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF][START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF]). The critical time T ⋆ is related to the measure of (0, 1) \ ω. Moreover, as a consequence of the Hilbert Uniqueness Method of J.-L. Lions [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF], the null controllability of ( 1) is equivalent to an observability inequality for the associated adjoint problem : there exists C > 0 such that

ϕ(•, 0), ϕ t (•, 0) 2 H ≤ C ϕ 2 L 2 (ω×(0,T )) , ∀(ϕ 0 , ϕ 1 ) ∈ H := L 2 (0, 1) × H -1 (0, 1) (4)
where (ϕ, ϕ 0 , ϕ 1 ) solves L ϕ = 0 in Q T , ϕ = 0 on Σ T ; (ϕ(•, 0), ϕ t (•, 0)) = (ϕ 0 , ϕ 1 ) in (0, 1).

(

) 5
We investigate in this work some questions related to the controllability of (1) for more general subsets q T ⊂ Q T where the support of the control function v depends on the time variable. A particular example is to consider the region between two curves, i.e. q T = (x, t) ∈ Q T ; a(t) < x < b(t), t ∈ (0, T ) .

A geometrical description is given by Figure 1.

To our knowledge, the control of PDEs with non-cylindrical support has been much less addressed in the literature. For the wave equation, we mention the contribution of Khapalov [START_REF] Khapalov | Controllability of the wave equation with moving point control[END_REF] providing observability results for a moving point sensor in the interior of the domain, allowing the author to avoid the usual difficulties related to strategic or non-strategic points. In particular, in the 1D setting, for any T > 0, the existence of controls continuous almost everywhere in (0, T ), supported over curves continuous almost everywhere is obtained for data in H 2 (0, 1) ∩ H 1 0 (0, 1) × H 1 0 (0, 1). More recently, let us mention two works concerning again the 1D wave equation both for initial data in H and any T > 2: the first one [START_REF] Castro | Exact controllability of the 1-D wave equation from a moving interior point[END_REF] analyzes the exact controllability from a moving interior point. By the way of the d'Alembert formula, an observability inequality is proved for a precise sets of curves {(γ(t), t)} t∈(0,T ) leading to moving controls in H -1 (∪ t∈(0,T ) γ(t) × {t}). The second one [START_REF] Cui | Exact controllability for a one-dimensional wave equation in non-cylindrical domains[END_REF] considers the controllability from the moving boundary of the form 1 + kt with k ∈ (0, 1], t ∈ (0, T ). In the case k < 1, the controllability in L 2 (0, T ) is proved by the way of the multiplier method: as is usual, a change of variable reduces the problem on fixed domains for wave equation with non constants coefficient. In the case k = 1 for which the speed of the moving endpoint is equal to the characteristic speed, the d'Alembert formula allows to characterize the reachable set. Eventually, we mention that in [START_REF] Martin | Null controllability of the structurally damped wave equation with moving control[END_REF], the authors prove that a mobile control support in time allows to recover the controllability of the damped wave equation y tty xxεy txx = 0 defined on the 1D torus, uniformly with respect to the parameter ε > 0 (due the appearance of an essential spectrum as ε → 0, the uniform controllability property can not be obtained if the control simply acts on a cylindrical domain).

Let us denote by C(a, b, T ) the class of domains of the form given by [START_REF] Castro | Exact controllability of the 1-D wave equation from a moving interior point[END_REF] for which the controllability holds, or equivalently the set of triplet (a, b, T ) ∈ C([0, T ], ]0, 1[) × C([0, T ], ]0, 1[) × R + for which the controllability hold. Obviously, this set is not empty: it suffices that T be large enough and that the domain {(x, t) ∈ Q T ; a(t) < x < b(t), t ∈ (0, T )} contains any rectangular domain (a 1 , b 1 ) × (0, T ) (i.e. that there exists a 1 > 0, b 1 > 0 such that a(t) ≤ a 1 < b 1 ≤ b(t) for all time t ∈ [0, T ]) and then to apply [START_REF] Bramble | Least-squares methods for Stokes equations based on a discrete minus one inner product[END_REF].

The first contribution of this work is the extension of the well-known observability inequality (4) to time-depending subsets q T , as for instance those given by [START_REF] Castro | Exact controllability of the 1-D wave equation from a moving interior point[END_REF]. Precisely, we prove that if T > 2 and any characteristic line starting at any point x ∈ (0, 1) at time t = 0 and following the optical geometric conditions when reflecting at the boundary x = 0, 1, meets the observation subset q T then the following estimate turns out to be true :

ϕ(•, 0), ϕ t (•, 0)) 2
H ≤ C ϕ 2 L 2 (q T ) + Lϕ 2 L 2 (0,T ;H -1 (0,1)) , ∀ϕ ∈ Φ [START_REF] Chapelle | The inf-sup test[END_REF] where Φ = {ϕ : ϕ ∈ L 2 (q T ), ϕ = 0 on Σ T such that Lϕ ∈ L 2 (0, T ; H -1 (0, 1))}. This so-called generalized inequality observability implies by duality arguments the null controllability of [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF] with controls in L 2 (q T ) for q T given by [START_REF] Castro | Exact controllability of the 1-D wave equation from a moving interior point[END_REF].

The second contribution of this work is the introduction of a numerical method for the approximation of the control of minimal L 2 (q T )-norm. Usually (for instance see [START_REF] Glowinski | Exact and approximate controllability for distributed parameter systems[END_REF][START_REF] Glowinski | Exact and approximate controllability for distributed parameter systems[END_REF][START_REF] Münch | A uniformly controllable and implicit scheme for the 1-D wave equation[END_REF]), such approximation is based on the minimization of the so-called conjugate functional :

min (ϕ0,ϕ1)∈H J ⋆ (ϕ 0 , ϕ 1 ) = 1 2 q T |ϕ| 2 dx dt+ < ϕ 1 , y 0 > H -1 (0,1),H 1 0 (0,1) - 1 0 ϕ 0 y 1 dx (8)
where (ϕ, ϕ 0 , ϕ 1 ) solves [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. Here and in the sequel, we use the following duality pairing :

< ϕ 1 , y 0 > H -1 (0,1),H 1 0 (0,1) = 1 0 ∂ x ((-∆) -1 ϕ 1 )(x) ∂ x y 0 (x) dx
where -∆ is the Dirichlet Laplacian in (0, 1). The well-posedness of this extremal problemprecisely the coercivity of J ⋆ over H -is given by the estimate [START_REF] Chapelle | The inf-sup test[END_REF]. The control of minimal L 2 (q T )norm is then given by v = -ϕ 1 q T where ϕ solves (5) associated to (ϕ 0 , ϕ 1 ), unique minimizer of J ⋆ . Since the domain q T evolves in time, the resolution of the wave equation by any usual time marching method leads to some technicalities (at the level of the numerical implementation). We therefore use, in the spirit of the works [START_REF] Cîndea | Numerical controllability of the wave equation through primal methods and carleman estimates[END_REF][START_REF] Cîndea | A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations[END_REF], an equivalent formulation where the time and the space variables are embedded in a time-space variational formulation, very appropriate to our non-cylindrical situation. Specifically, we consider the following extremal problem over the variable ϕ:

min ϕ∈W Ĵ⋆ (ϕ) = 1 2 q T |ϕ| 2 dx dt+ < ϕ 1 , y 0 > H -1 (0,1),H 1 0 (0,1) - 1 0 ϕ 0 y 1 dx, (9)
with W = ϕ : ϕ ∈ L 2 (q T ), ϕ = 0 on Σ T such that Lϕ = 0 ∈ L 2 (0, T ; H -1 (0, 1)) . The unknown is now the whole solution ϕ, constrained by the equality Lϕ = 0 in L 2 (0, T ; H -1 (0, 1)). This constraint is taken into account with a Lagrangian multiplier λ ∈ L 2 (0, T ; H 1 0 (0, 1)) through a mixed formulation in (ϕ, λ) solved using a conformal finite element approximation in time and space.

This paper is organized as follows. In Section 2, we prove the generalized observability inequality (7) over the Hilbert space Φ (see Proposition 2.1) leading to the uniform controllability of (1) in V . In Section 3, we adapt [START_REF] Cîndea | A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations[END_REF] and introduce the mixed formulation [START_REF] Münch | Optimal distribution of the internal null control for the onedimensional heat equation[END_REF] in the variable (ϕ, λ) ∈ Φ × L 2 (0, T ; H 1 0 (0, 1)) equivalent to the extremal problem [START_REF] Cîndea | Numerical controllability of the wave equation through primal methods and carleman estimates[END_REF]. The variable λ can be interpreted as a Lagrange multiplier for the equality constraint Lϕ = 0 in L 2 (0, T ; H -1 (0, 1)). We employ the estimate [START_REF] Chapelle | The inf-sup test[END_REF] to prove the well-posedness of this mixed formulation. In particular, we prove an inf-sup condition for the pair (ϕ, λ). Moreover, it turns out that the multiplier λ, unique solution of the mixed formulation, coincides with the controlled state y, solution of (1) (in the weak-sense) (see section 3.1). This property allows to define Section 3.2 an another equivalent extremal problem (the so-called primal problem, dual of the problem ( 9)) in the controlled solution y only (see Proposition 3.1), without the introduction of any penalty parameter. The corresponding elliptic problem in L 2 (0, T ; H 1 0 (0, 1)) is solved by the way of a conjugate gradient algorithm. Section 4 is devoted to the numerical approximation of the mixed formulation as well as some numerical experiments. We emphasize the robustness of the approach leading notably to the strong convergence of discrete sequence {v h } toward the controls for various geometries of q T . Eventually, Section 5 concludes with some perspectives: in particular, we highlight the natural extension of this work which consists in optimizing the control of (1) with respect to the support q T . 0

1 2 t 0 0.5 1 x γ(t) δ 0 q T Q T 0 1 2 t 0 0.5 1 x γ(t) q T Q T Figure 1: Time dependent domains q T included in Q T .

A generalized observability inequality

Let q T ⊂ (0, 1) × (0, T ) be an open set. We define the vectorial space Φ = ϕ : ϕ ∈ L 2 (q T ), ϕ = 0 on Σ T such that Lϕ ∈ L 2 (0, T ; H -1 (0, 1)) .

Endowed with the following inner product

(ϕ, ϕ) Φ = q T ϕ(x, t)ϕ(x, t) dx dt + η T 0 < Lϕ, Lϕ > H -1 (0,1),H -1 (0,1) dt, (10)
for any fixed η > 0, the space Φ is an Hilbert space.

In this section, we prove the following result.

Proposition 2.1 Assume that T > 2 and q T ⊂ (0, 1) × (0, T ) is a finite union of connected open sets and satisfies the following hypotheses: Any characteristic line starting at the point x ∈ (0, 1) at time t = 0 and following the optical geometric laws when reflecting at the boundaries x = 0, 1 must meet q T . Then, there exists C > 0 such that the following estimate holds :

ϕ(•, 0), ϕ t (•, 0)) 2 H ≤ C ϕ 2 L 2 (q T ) + Lϕ 2 L 2 (0,T ;H -1 (0,1)) , ∀ϕ ∈ Φ. ( 11
)
Proof: We proceed in several steps:

Step 1: Let ϕ be a smooth solutions of the wave equation [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. Note that ϕ can be extended in a unique way to a function, still denoted by ϕ, in (x, t) ∈ (0, 1) × R satisfying Lϕ = 0 and the boundary conditions ϕ(0, t) = ϕ(1, t) = 0 for all t ∈ R. In this first step we prove that for such extension the following holds: For each t ∈ R, x ∈ (0, 1) and δ > 0

t+δ t-δ ϕ x (0, s) ds ≤ 1 δ U δ (x,t+x) |ϕ x (y, s)| 2 + |ϕ t (y, s)| 2 dy ds, ( 12
) t+δ t-δ ϕ x (0, s) ds ≤ 1 δ U δ (x,t-x) |ϕ x (y, s)| 2 + |ϕ t (y, s)| 2 dy ds, (13)
where U δ (x,t) is a neighborhood of (x, t) of the form U δ (x,t) = {(y, s) such that |x -y| + |t -s| < δ}. Let us start with [START_REF] Daniel | The approximate minimization of functionals[END_REF]. Note that wave equation is symmetric with respect to the time and space variables. Therefore, the D'Alembert formulae can be used changing the time and space role, i.e.

1 2

t+x t-x u x (0, s) ds = u(x, t), (x, t) ∈ (0, 1) × R, (14)
where we have taken into account the boundary condition u(0, t) = 0. Consider now x = x 0t in ( 14) and differentiate with respect to time. Then,

-ϕ x (0, 2t -x 0 ) = -ϕ x (x 0 -t, t) + ϕ t (x 0 -t, t),
that written in the original variables (x, t) gives,

ϕ x (0, t -x) = ϕ x (x, t) -ϕ t (x, t), or equivalently ϕ x (0, t) = ϕ x (x, t + x) -ϕ t (x, t + x), t ∈ R, x ∈ (0, 1). ( 15
)
Integrating the square of (15) in (y, s) ∈ U δ (x,t+x) with the parametrization

y = x + u-v √ 2 s = t + x + u+v √ 2 , |u|, |v| < δ/ √ 2, we obtain δ/ √ 2 -δ/ √ 2 δ/ √ 2 -δ/ √ 2 |u x (0, t + 2v/ √ 2)| 2 dudv = U δ x,t+x |u x (y, s) -u t (y, s)| 2 dyds.
Therefore, with the change s = t + 2v/ √ 2 in the first integral, we easily obtain [START_REF] Daniel | The approximate minimization of functionals[END_REF]. Formula ( 13) is obtained in a similar way but this time with the change x = x 0 + t in [START_REF] Fortin | Augmented Lagrangian methods[END_REF].

Step 2. Here we show that there exists a constant C > 0 such that,

ϕ(•, 0), ϕ t (•, 0)) 2 V ≤ C ϕ t 2 L 2 (q T ) + ϕ x 2 L 2 (q T ) , (16)
for any ϕ ∈ W and initial data in V . We may assume that ϕ is smooth since the general case can be obtained by a usual density argument. We also assume that ϕ is extended to (x, t) ∈ (0, 1) × R by assuming that it satisfies the wave equation and boundary conditions in this region. This extension is unique and 2-periodic in time. The region q T is also extended to qT to take advantage of the time periodicity of the solution ϕ. We define, qT = k∈Z Z {(x, t) such that (x, t + 2k) ∈ q T } .

The key point now is to observe that the hypotheses on q T , namely the fact that any characteristic line starting at point (x, 0) and following the optical geometric laws when reflecting at the boundary must meet q T , is equivalent to the following: For any point (0, t) with t ∈ [0, T ] there exists one characteristic line (either (x, t + x) with x ∈ (0, 1) or (x, tx)) that meets qT .

Thus, given t ∈ [0, T ] we can apply either [START_REF] Daniel | The approximate minimization of functionals[END_REF] with (x, t + x) ∈ qT or [START_REF] Dunavant | High degree efficient symmetrical Gaussian quadrature rules for the triangle[END_REF] with (x, tx) ∈ qT . Moreover, as qT is an open set, we may assume δ sufficiently small so that the neighborhoods U δ (x,t+x) ⊂ qT or U δ (x,t-x) ⊂ qT . In particular we see that for any t ∈ [0, T ] there exists δ t > 0 and

C t > 0 such that t+δt t-δt |ϕ x (0, s)| 2 ds ≤ C t qT |ϕ t | 2 + |ϕ x | 2 dxdt.
By compacity, there exists a finite number of times t 1 , ..., t n such that ∪ i=1,..,n (t iδ ti , t i + δ ti ) covers the whole interval [0, T ] and therefore, by adding the corresponding inequalities, there exists C > 0 such that

T 0 |ϕ x (0, s)| 2 ds ≤ C qT |ϕ t | 2 + |ϕ x | 2 dxdt. ( 17
)
The fact that we can replace qT by q T is easily checked by the 2-periodicity of ϕ in time. Finally, it remains to see that we can replace the left hand side in [START_REF] Glowinski | Exact and approximate controllability for distributed parameter systems[END_REF] by ϕ(•, 0), ϕ t (•, 0)) 2 V , but this is a consequence of the well-known boundary observability inequality (see for instance [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF])

ϕ(•, 0), ϕ t (•, 0)) 2 V ≤ C T 0 |ϕ x (0, s)| 2 ds
which holds for T > 2.

Step 3. Here we show that we can substitute ϕ x by ϕ in the right hand side of ( 16), i.e.

ϕ(•, 0), ϕ t (•, 0)) 2 V ≤ C ϕ t 2 L 2 (q T ) + ϕ 2 L 2 (q T ) , (18)
for any ϕ ∈ W and initial data in V . In fact, this requires to extend slightly the observation zone q T . Instead, we observe that if q T satisfies the hypotheses in Proposition 2.1 then there exists a smaller open subset qT ⊂ q T that still satisfies the same hypotheses and such that the closure of qT is included in q T . Thus, (16) must hold as well for qT . Let us introduce now a function η ≥ 0 which satisfies the following hypotheses:

η ∈ C 1 ((0, 1) × (0, T )), supp(η) ⊂ q T , η t L ∞ + η 2 x /η L ∞ ≤ C 1 in q T η > η 0 > 0 in qT , with η 0 > 0 constant.
As q T is a finite union of connected open sets, the function η can be easily obtained by convolution of the characteristic function of qT with a positive mollifier.

Multiplying the equation of ϕ by ηϕ and integrating by parts we easily obtain

q T η|ϕ x | 2 dx dt = q T η|ϕ t | 2 dx dt + q T (η t ϕϕ t -η x ϕϕ x ) dx dt ≤ q T η|ϕ t | 2 dxdt + η t L ∞ (q T ) 2 q T (|ϕ| 2 + |ϕ t |) dx dt + 1 2 q T ( η 2 x η ϕ 2 + ηϕ 2 x ) dx dt.
Therefore,

q T η|ϕ x | 2 dx dt ≤ C q T (|ϕ t | 2 + |ϕ| 2 ) dx dt,
for some constant C > 0, and we obtain

ϕ x 2 L 2 (q T ) ≤ C -1 2 q T η|ϕ x | 2 dx dt ≤ C -1 2 C q T (|ϕ t | 2 + |ϕ| 2 ) dx dt.
This combined with [START_REF] Glowinski | Exact and approximate controllability for distributed parameter systems[END_REF] for qT provides [START_REF] Khapalov | Controllability of the wave equation with moving point control[END_REF].

Step 4. Here we prove that we can remove the second term in the right hand side of (18), i.e.

ϕ(•, 0), ϕ t (•, 0)) 2 V ≤ C ϕ t 2 L 2 (q T ) , (19)
for any ϕ ∈ W and initial data in V . Note that, for each time t ∈ [0, T ] and each ω ⊂ Ω we have the following regularity estimate

b(t) a(t) |ϕ(x, t)| 2 dx ≤ ϕ(•, 0), ϕ t (•, 0)) 2 H , for all t ∈ [0, T ]
Therefore, integrating in time, we easily obtain

ϕ 2 L 2 (q T ) ≤ T ϕ(•, 0), ϕ t (•, 0)) 2 H .
We now substitute this inequality in ( 16)

ϕ(•, 0), ϕ t (•, 0)) 2 V ≤ C ϕ t 2 L 2 (q T ) + ϕ(•, 0), ϕ t (•, 0)) 2 H . (20)
Inequality ( 19) is finally obtained by contradiction. Assume that it is not true. Then, there exists a sequence (ϕ k (•, 0),

ϕ k t (•, 0))) k>0 ∈ V such that ϕ k (•, 0), ϕ k t (•, 0)) 2 V = 1, ∀k > 0 ϕ k t 2 L 2 (q T ) → 0, as k → ∞.
Therefore, there exists a subsequence, still denoted by the index k, such that (ϕ k (•, 0),

ϕ k t (•, 0)) → (ϕ ⋆ (•, 0), ϕ ⋆ t (•, 0)
) weakly in V and strongly in H (by the compactness of the inclusion V ⊂ H). Passing to the limit in the equation we easily see that the solution associated to (ϕ ⋆ (•, 0), ϕ ⋆ t (•, 0)), ϕ ⋆ must vanish at q T and therefore, by [START_REF] Khapalov | Controllability of the wave equation with moving point control[END_REF], ϕ ⋆ = 0. In particular this is in contradiction with [START_REF] Martin | Null controllability of the structurally damped wave equation with moving control[END_REF] since it implies that the right hand side in [START_REF] Martin | Null controllability of the structurally damped wave equation with moving control[END_REF] vanishes as k → ∞ while the left hand side is bounded by below.

Step 5. We now write [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF] with respect to the weaker norm. In particular, we obtain

ϕ(•, 0), ϕ t (•, 0)) 2 H ≤ C ϕ 2 L 2 (q T ) , (21)
for any ϕ ∈ Φ with Lϕ = 0. Let η ∈ Φ be the solution of Lη = 0 and initial data (∆ -1 ϕ t (•, 0), ϕ(•, 0)) ∈ V where ∆ designates the Dirichlet Laplacian in (0, 1). Let us write η(x, t) = η(x, 0) + t 0 ϕ(x, s) ds, for all (x, t) ∈ Q T . Then, inequality [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF] on η and the fact that ∆ is an isomorphism from H 1 0 (0, 1) to

H -1 (0, 1), provide (ϕ(•, 0), ϕ t (•, 0), ) 2 H = (∆ -1 ϕ t (•, 0), ϕ(•, 0)) 2 V ≤ C η t 2 L 2 (q T ) = C ϕ 2 L 2 (q T ) .
Step 6. Here we finally obtain [START_REF] Cui | Exact controllability for a one-dimensional wave equation in non-cylindrical domains[END_REF]. Given ϕ ∈ Φ we can decompose it as ϕ = ϕ 1 + ϕ 2 where ϕ 1 , ϕ 2 ∈ Φ solve

Lϕ 1 = Lϕ, ϕ 1 (•, 0) = (ϕ 1 ) t (•, 0) = 0 Lϕ 2 = 0, ϕ 2 (•, 0) = ϕ(•, 0), (ϕ 2 ) t (•, 0) = ϕ t (•, 0).
From Duhamel's principle, we can write

ϕ 1 (•, t) = t 0 ψ(•, t -s, s)ds
where ψ(x, t, s) solves, for each value of the parameter s ∈ (0, t),

Lψ(•, •, s) = 0, ψ(•, 0, s) = 0, ψ t (•, 0, s) = Lϕ(•, s).
Therefore,

ϕ 1 2 L 2 (q T ) ≤ T 0 ψ(•, •, s) 2 L 2 (q T ) ds ≤ C T 0 ψ(•, 0, s), ψ t (•, 0, s)) 2 H ds ≤ C Lϕ 2 L 2 (0,T ;H -1 (0,1)) (22)
Combining [START_REF] Münch | A uniformly controllable and implicit scheme for the 1-D wave equation[END_REF] and estimate ( 21) for ϕ 2 we obtain

ϕ(•, 0), ϕ t (•, 0)) 2 H = ϕ 2 (•, 0), (ϕ 2 ) t (•, 0)) 2 H ≤ C ϕ 2 2 L 2 (q T ) ≤ C ϕ 2 L 2 (q T ) + ϕ 1 2 L 2 (q T ) ≤ C ϕ 2 L 2 (q T ) + Lϕ 2 L 2 (0,T ;H -1 ) .
✷ Remark 1 1. The hypotheses on q T stated in Proposition 2.1 are optimal in the following sense: If there exists a subinterval ω 0 ⊂ (0, 1) for which all characteristics starting in ω 0 and following the geometrical optics conditions when getting to the boundary x = 0, 1, do not meet q T , then the inequality fails to hold. This is easily seen by considering particular solutions of the wave equation which initial data supported in ω 0 .

2. The proof of inequality [START_REF] Cui | Exact controllability for a one-dimensional wave equation in non-cylindrical domains[END_REF] above does not provide an estimate on the dependence of the constant with respect to q T .

3. In the cylindrical situation, i.e. q T = (a, b) × (0, T ), a generalized Carleman inequality, valid for the wave equation with variable coefficients, have been obtained in [START_REF] Cîndea | Numerical controllability of the wave equation through primal methods and carleman estimates[END_REF] (see also [START_REF] Yao | On the observability inequalities for exact controllability of wave equations with variable coefficients[END_REF]). To our knowledge, the extension of Proposition 2.1 to the wave equation with variable coefficients is still open and a priori can not be obtained by the method used in this section.

3 Control of minimal L 2 (q T )-norm: a mixed reformulation

We now adapt in this section the work [START_REF] Cîndea | A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations[END_REF] and present a mixed formulation based on the optimality conditions associated to the extremal problem (8) (section 3.1). From a numerical point of view, this mixed formulation is more appropriate to the non-cylindrical situation considered in this work.

Mixed reformulation of the controllability problem

As described at length in [START_REF] Cîndea | A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations[END_REF], the starting point of the reformulation is the dual problem [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. Thus, in order to avoid the minimization of the functional J ⋆ with respect to the initial data (ϕ 0 , ϕ 1 ), we now present a direct way to approximate the control of minimal square integrable norm, in the spirit of the primal approach developed in [START_REF] Cîndea | Numerical controllability of the wave equation through primal methods and carleman estimates[END_REF]. Since the variable ϕ, solution of (5), is completely and uniquely determined by the data (ϕ 0 , ϕ 1 ), the main idea of the reformulation is to keep ϕ as main variable and consider the following extremal problem:

min ϕ∈W Ĵ⋆ (ϕ) = 1 2 q T |ϕ| 2 dx dt+ < ϕ t (•, 0), y 0 > H -1 (0,1),H 1 0 (0,1) - 1 0 ϕ(•, 0) y 1 dx, (23)
where W = ϕ : ϕ ∈ L 2 (q T ), ϕ = 0 on Σ T such that Lϕ = 0 , endowed with the same inner product than Φ, given in [START_REF] Cîndea | A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations[END_REF].

Remark that from [START_REF] Cui | Exact controllability for a one-dimensional wave equation in non-cylindrical domains[END_REF] the property ϕ ∈ W implies that (ϕ(•, 0), ϕ t (•, 0)) ∈ H, so that the functional Ĵ⋆ is well-defined over W . Therefore, the minimization of Ĵ⋆ is evidently equivalent to the minimization of J ⋆ over H.

The main variable is now ϕ submitted to the constraint equality Lϕ = 0 as an L 2 (0, T ; H -1 (0, 1)) function. This constraint is addressed introducing a Lagrangian multiplier λ ∈ L 2 (0, T ; H 1 0 (Ω)) as follows:

We consider the following problem : find (ϕ, λ)

∈ Φ × L 2 (0, T ; H 1 0 (0, 1)) solution of a(ϕ, ϕ) + b(ϕ, λ) = l(ϕ), ∀ϕ ∈ Φ b(ϕ, λ) = 0, ∀λ ∈ L 2 (0, T ; H 1 0 (0, 1)), (24)
where

a : Φ × Φ → R, a(ϕ, ϕ) = q T ϕ ϕ dx dt b : Φ × L 2 (0, T ; H 1 0 (0, 1)) → R, b(ϕ, λ) = T 0 < Lϕ, λ > H -1 (0,1),H 1 0 (0,1) dt = Q T ∂ x (-∆ -1 (Lϕ)) • ∂ x λ dx dt l : Φ → R, l(ϕ) = -< ϕ t (•, 0), y 0 > H -1 (0,1),H 1 0 (0,1) + 1 0 ϕ(•, 0) y 1 dx.
We have the following result :

Theorem 3.1 Under the hypothesis of Proposition 2.1, we have :

(i) The mixed formulation (24) is well-posed.

(ii) The unique solution (ϕ, λ)

∈ Φ × L 2 (0, T ; H 1 0 (0, 1)) is the unique saddle-point of the La- grangian L : Φ × L 2 (0, T ; H 1 0 (0, 1)) → R defined by L(ϕ, λ) = 1 2 a(ϕ, ϕ) + b(ϕ, λ) -l(ϕ). (25)
(iii) The optimal function ϕ is the minimizer of Ĵ⋆ over Φ while the optimal function λ ∈ L 2 (0, T ; H 1 0 (0, 1)) is the state of the controlled wave equation (1) in the weak sense (associated to the control -ϕ 1 q T ).

Proof -We easily check that the bilinear form a is continuous over Φ × Φ, symmetric and positive and that the bilinear form b is continuous over Φ × L 2 (0, T ; H 1 0 (0, 1)). Furthermore, the continuity of the linear form l over Φ is a direct consequence of the generalized observability inequality [START_REF] Cui | Exact controllability for a one-dimensional wave equation in non-cylindrical domains[END_REF]:

|l(ϕ)| ≤ (y 0 , y 1 ) V C max(1, η -1 ) ϕ Φ , ∀ϕ ∈ Φ.
Therefore, the well-posedness of the mixed formulation is a consequence of the following two properties (see [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]):

• a is coercive on N (b), where N (b) denotes the kernel of b :

N (b) = {ϕ ∈ Φ such that b(ϕ, λ) = 0 for every λ ∈ L 2 (0, T ; H 1 0 (0, 1))};
• b satisfies the usual "inf-sup" condition over Φ × L 2 (0, T ; H 1 0 (0, 1)): there exists δ > 0 such that inf λ∈L 2 (0,T ;H 1 0 (0,1))

sup ϕ∈Φ b(ϕ, λ) ϕ Φ λ L 2 (0,T ;H 1 0 (0,1)) ≥ δ. ( 26
)
From the definition of a, the first point is clear : for all ϕ ∈ N (b) = W , a(ϕ, ϕ) = ϕ 2 W . Let us check the inf-sup condition [START_REF] Yao | On the observability inequalities for exact controllability of wave equations with variable coefficients[END_REF]. For any fixed λ 0 ∈ L 2 (0, T ; H 1 0 (0, 1)), we define the (unique) element ϕ 0 such that Lϕ 0 = -∆λ 0 in Q T and such that ϕ 0 (•, 0) = 0 in L 2 (0, 1) and ϕ 0,t (•, 0) = 0 in H -1 (0, 1). ϕ 0 is therefore solution of the wave equation with source term -∆λ 0 ∈ L 2 (0, T ; H -1 (0, 1)), null Dirichlet boundary condition and zero initial state. We then use the following estimate (see for instance Chapter 1 in [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF]): there exists a constant C Ω,T > 0 such that

ϕ 0 L 2 (Q T ) ≤ C Ω,T -∆λ 0 L 2 (0,T ;H -1 (0,1)) ≤ C Ω,T λ 0 L 2 (0,T ;H 1 0 (0,1)) . (27)
Consequently, ϕ 0 ∈ Φ. In particular, we have b(ϕ

0 , λ 0 ) = λ 0 2 L 2 (0,T ;H 1 0 (0,1)) and sup ϕ∈Φ b(ϕ, λ 0 ) ϕ Φ λ 0 L 2 (Q T ) ≥ b(ϕ 0 , λ 0 ) ϕ 0 Φ λ 0 L 2 (Q T ) = λ 0 2 L 2 (0,T ;H 1 0 (0,1)) ϕ 0 2 L 2 (q T ) + η λ 0 2 L 2 (0,T ;H 1 0 (0,1)) 1 2 λ 0 L 2 (0,T ;H 1 0 (0, 1))
.

Combining the above two inequalities, we obtain

sup ϕ0∈Φ b(ϕ 0 , λ 0 ) ϕ 0 Φ λ 0 L 2 (0,T ;H 1 0 (0,1)) ≥ 1 C 2 Ω,T + η
and, hence, [START_REF] Yao | On the observability inequalities for exact controllability of wave equations with variable coefficients[END_REF] holds with δ = C 2 Ω,T + η

-1
2 . The point (ii) is due to the symmetry and to the positivity of the bilinear form a. (iii). The equality b(ϕ, λ) = 0 for all λ ∈ L 2 (0, T ; H 1 0 (0, 1)) implies that Lϕ = 0 as an L 2 (0, T ; H -1 (0, 1)) function, so that if (ϕ, λ) ∈ Φ × L 2 (0, T ; H 1 0 (0, 1)) solves the mixed formulation, then ϕ ∈ W and L(ϕ, λ) = Ĵ⋆ (ϕ). Finally, the first equation of the mixed formulation reads as follows :

q T ϕ ϕ dx dt + T 0 < Lϕ, λ > H -1 ,H 1 0 dt = l(ϕ), ∀ϕ ∈ Φ, or equivalently, since the control of minimal L 2 (q T ) norm is given by v = -ϕ 1 q T , Q T -v1 q T ϕ dx dt + T 0 < Lϕ, λ > H -1 ,H 1 0 dt = l(ϕ), ∀ϕ ∈ Φ.
But this means that λ ∈ L 2 (0, T, H 1 0 (0, 1)) is solution of the wave equation in the transposition sense. Since (y 0 , y 1 ) ∈ V and v ∈ L 2 (q T ), λ must coincide with the unique weak solution to [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF]. ✷ Therefore, Theorem 3.1 reduces the search of the control of square minimal norm to the resolution of the mixed formulation [START_REF] Münch | Optimal distribution of the internal null control for the onedimensional heat equation[END_REF], or equivalently to the search of the saddle point for L. In general, it is very convenient (and actually necessary in Section 3.2) to "augment" the Lagrangian (see [START_REF] Fortin | Augmented Lagrangian methods[END_REF]), and consider instead the Lagrangian L r defined for any r > 0 by

   L r (ϕ, λ) := 1 2 a r (ϕ, ϕ) + b(ϕ, λ) -l(ϕ), a r (ϕ, ϕ) := a(ϕ, ϕ) + r L ϕ 2 L 2 (0,T ;H -1 (0,1)) .
Since a(ϕ, ϕ) = a r (ϕ, ϕ) on W , the Lagrangian L and L r share the same saddle-point.

Remark 2

The estimate [START_REF] Cui | Exact controllability for a one-dimensional wave equation in non-cylindrical domains[END_REF] may also be used to extend the work [START_REF] Cîndea | Numerical controllability of the wave equation through primal methods and carleman estimates[END_REF] to the non-cylindrical situation. Reference [START_REF] Cîndea | Numerical controllability of the wave equation through primal methods and carleman estimates[END_REF] considers the pair (y, v) solution of (1-3) which minimize the following L 2 -weighted functional

J(y, v) := 1 2 Q T ρ 2 (x, t)|y| 2 dx dt + 1 2 q T ρ 2 0 (x, t)|v| 2 dx dt for any weights ρ, ρ 0 ∈ C(Q T , R ⋆ + ).
Assuming |ρ| and |ρ 0 | uniformly positive by below, the unique minimizer (y, v) is expressed in term of the auxiliary variable p ∈ P := {p :

ρ -1 Lp ∈ L 2 (Q T ), ρ -1 0 p ∈ L 2 (q T ), p = 0 on Σ T } as follows : y = -ρ -2 Lp, v = ρ -2 0 p 1 q T on Q T
where p is the solution of the variational formulation

Q T ρ -2 Lp Lq dx dt+ q T ρ -2 0 p q dx dt = 1 0 y 1 q(•, 0)dx-< y 0 , q t (•, 0) > H 1 0 (0,1),H -1 (0,1) , ∀q ∈ P.
The well-posedness of this formulation is given by the estimate (11).

Dual problem of the extremal problem (23)

The mixed formulation allows to solve simultaneously the dual variable ϕ, argument of the conjugate functional [START_REF]Optimal design of the support of the control for the 2-D wave equation: a numerical method[END_REF], and the Lagrangian multiplier λ. Since λ turns out to be the controlled state of (1), we may qualify λ as the primal variable of the controllability problem. We derive in this section the corresponding extremal problem involving only that variable λ. For any r > 0, let us define the linear operator A r from L 2 (0, T ; H 1 0 (0, 1)) into L 2 (0, T ; H 1 0 (0, 1)) by A r λ := -∆ -1 (Lϕ), ∀λ ∈ L 2 (0, T ; H 1 0 (0, 1)) where ϕ ∈ Φ is the unique solution to

a r (ϕ, ϕ) = b(ϕ, λ), ∀ϕ ∈ Φ. ( 28
)
Notice that the assumption r > 0 is needed here in order to guarantee the well-posedness of (28). Precisely, for any r > 0, the form a r defines a norm equivalent to the norm on ϕ.

We have the following important lemma :

Lemma 3.1 For any r > 0, the operator A r is a strongly elliptic, symmetric isomorphism from

L 2 (0, T ; H 1 0 (0, 1)) into L 2 (0, T ; H 1 0 (0, 1)).
Proof-From the definition of a r , we easily get that A r λ L 2 (0,T ;

H 1 0 (0,1)) ≤ r -1 λ L 2 (0,T ;H 1 0 (0,1))
and the continuity of A r . Next, consider any λ ′ ∈ L 2 (0, T ; H 1 0 (0, 1)) and denote by ϕ ′ the corresponding unique solution of (28) so that A r λ ′ := -∆ -1 (Lϕ ′ ). Relation (28) with ϕ = ϕ ′ then implies that

T 0 < A r λ ′ , λ > H 1 0 (0,1),H 1 0 (0,1) dt = a r (ϕ, ϕ ′ ) ( 29
)
and therefore the symmetry and positivity of A r . The last relation with λ ′ = λ and the estimate [START_REF] Cui | Exact controllability for a one-dimensional wave equation in non-cylindrical domains[END_REF] imply that A r is also positive definite.

Finally, let us check the strong ellipticity of A r , equivalently that the bilinear functional (λ, λ ′ ) → T 0 < A r λ, λ ′ > H 1 0 (0,1),H 1 0 (0,1) dt is L 2 (0, T ; H 1 0 (0, 1))-elliptic. Thus we want to show that

T 0 < A r λ, λ > H 1 0 (0,1),H 1 0 (0,1) dt ≥ C λ 2 L 2 (0,T ;H 1 0 (0,1)) , ∀λ ∈ L 2 (0, T ; H 1 0 (0, 1)) ( 30
)
for some positive constant C. Suppose that (30) does not hold; there exists then a sequence {λ n } n≥0 of L 2 (0, T ; H 1 0 (0, 1)) such that λ n L 2 (0,T,H 1 0 (0,1)) = 1, ∀n ≥ 0, and lim

n→∞ T 0 < A r λ n , λ n > H 1 0 (0,1),H 1 0 (0,1) dt = 0.
Let us denote by ϕ n the solution of (28) corresponding to λ n . From (29), we then obtain that lim n→∞ Lϕ n L 2 (0,T,H -1 (0,1)) = 0, lim

n→+∞ ϕ n L 2 (q T ) = 0 (31)
and thus lim n→∞

T 0 < -∆ -1 (Lϕ), λ n > H 1 0 (0,1),H 1
0 (0,1) dt = 0 for all ϕ ∈ Φ (and so the L 2 (0, T ; H 1 0 (0, 1))weak-convergence of λ n toward 0).

From (28) with ϕ = ϕ n and λ n , we have

T 0 < -r∆ -1 (Lϕ n ) -λ n , -∆ -1 (Lϕ) > H 1 0 (0,1),H 1 0 (0,1) dt + q T ϕ n ϕ dx dt = 0, ∀ϕ ∈ Φ. ( 32
)
We define the sequence {ϕ n } n≥0 as follows :

     Lϕ n = r Lϕ n + ∆λ n , in Q T , ϕ n (0, •) = ϕ n (1, •) = 0, in (0, T ), ϕ n (•, 0) = ϕ n,t (•, 0) = 0, in (0, 1)
so that, for all n, ϕ n is the solution of the wave equation with zero initial data and source term r Lϕ n + ∆λ n in L 2 (0, T ; H -1 (0, 1)). Using again (27), we get ϕ n L 2 (q T ) ≤ C Ω,T rLϕ n + ∆λ n L 2 (0,T ;H -1 (0,1)) , so that ϕ n ∈ Φ. Then, using (32), we get

-r∆ -1 (Lϕ n ) -λ n L 2 (0,T ;H 1 0 (0,1)) ≤ C Ω,T ϕ n L 2 (q T ) .
Then, from (31), we conclude that lim n→+∞ λ n L 2 (0,T ;H 1 0 (0,1)) = 0 leading to a contradiction and to the strong ellipticity of the operator A r .

✷ The introduction of the operator A r is motivated by the following proposition :

Proposition 3.1 Let ϕ 0 ∈ Φ the unique solution of a r (ϕ 0 , ϕ) = l(ϕ), ∀ϕ ∈ Φ and let J ⋆⋆ : L 2 (0, T ; H 1 0 (0, 1)) → L 2 (0, T ; H 1 0 (0, 1
)) the functional defined by

J ⋆⋆ (λ) = 1 2 T 0 < A r λ, λ > H 1 0 (0,1),H 1 0 (0,1) dt -b(ϕ 0 , λ).
The following equality holds :

sup λ∈L 2 (0,T ;H 1 0 (0,1)) inf ϕ∈Φ L r (ϕ, λ) = - inf λ∈L 2 (0,T ;H 1 0 (0,1)) J ⋆⋆ (λ) + L r (ϕ 0 , 0).
Proof-For any λ ∈ L 2 (0, T ; H 1 0 (0, 1)), let us denote by ϕ λ ∈ Φ the minimizer of ϕ → L r (ϕ, λ). ϕ λ satisfies the equation a r (ϕ λ , ϕ) + b(ϕ, λ) = l(ϕ), ∀ϕ ∈ Φ and can be decomposed as follows :

ϕ λ = ψ λ + ϕ 0 where ψ λ ∈ Φ solves a r (ψ λ , ϕ) + b(ϕ, λ) = 0, ∀ϕ ∈ Φ.
We then have

inf ϕ∈Φ L r (ϕ, λ) = L r (ϕ λ , λ) = L r (ψ λ + ϕ 0 , λ) = 1 2 a r (ψ λ + ϕ 0 , ψ λ + ϕ 0 ) + b(ψ λ + ϕ 0 , λ) -l(ψ λ + ϕ 0 ) := X 1 + X 2 + X 3 with      X 1 = 1 2 a r (ψ λ , ψ λ ) + b(ψ λ , λ) + b(ϕ 0 , λ) X 2 = a r (ψ λ , ϕ 0 ) -l(ψ λ ), X 3 = 1 2 a r (ϕ 0 , ϕ 0 ) -l(ϕ 0 ).
From the definition of ϕ 0 , X 2 = 0 while X 3 = L r (ϕ 0 , 0). Eventually, from the definition of ψ λ ,

X 1 = - 1 2 a r (ψ λ , ψ λ ) + b(ϕ 0 , λ) = - 1 2 T 0 < A r λ, λ > H 1 ,H 1 dt + b(ϕ 0 , λ)
and the result follows. ✷ From the ellipticity of the operator A r , the minimization of the functional J ⋆⋆ over L 2 (0, T, H 1 0 ) is well-posed. It is interesting to note that with this extremal problem involving only λ, we are coming to the primal variable, controlled solution of (1) (see Theorem 3.1, (iii)). Due to the constraint (3), the direct minimization of the null controllability problem by a penalty method with respect to the controlled state is usually avoided in practice. Here, any constraint equality is assigned to the variable λ.

From the symmetry and ellipticity of the operator A r , the conjugate gradient algorithm is very appropriate to minimize J ⋆⋆ , and consequently solve the mixed formulation [START_REF] Münch | Optimal distribution of the internal null control for the onedimensional heat equation[END_REF]. The conjugate gradient algorithm reads as follows :

(i) Let λ 0 ∈ L 2 (0, T ; H 1 0 (0, 1)) be a given function.

(ii) Compute ϕ 0 ∈ Φ solution to a r (ϕ 0 , ϕ) + b(ϕ, λ 0 ) = l(ϕ), ∀ϕ ∈ Φ and g 0 = -∆ -1 (L ϕ 0 ) then set w 0 = g 0 .

(iii) For n ≥ 0, assuming that λ n , g n and w n are known, compute ϕ n ∈ Φ solution to

a r (ϕ n , ϕ) = b(ϕ, w n ), ∀ϕ ∈ Φ and g n = -∆ -1 (Lϕ n ) and then ρ n = g n 2
L 2 (0,T ;H 1 0 (0,1)) /(g n , w n ) L 2 (0,T ;H 1 0 (0,1)) . Update λ n and g n by

λ n+1 = λ n -ρ n w n , g n+1 = g n -ρ n g n . If g n+1 L 2 (0,T ;H 1 0 (0,1)) / g 0 L 2 (0,T ;H 1 0 (0,1)) ≤ ε, take λ = λ n+1 . Else, compute γ n = g n+1 2 L 2 (0,T ;H 1 0 (0,1)) / g n 2 L 2 (0,T ;H 1 0 (0,1))
and update w n via w n+1 = g n+1 + γ n w n .

Do n = n + 1 and return to step (iii).

As mentioned in [START_REF] Glowinski | Handbook of numerical analysis[END_REF] where this approach is discussed at length for Navier-Stokes type systems, this algorithm can be viewed as a sophisticated version of Arrow-Hurwicz-Uzawa type method.

Concerning the speed of convergence of the conjugate gradient algorithm (i)-(iii), it follows from for instance [START_REF] Daniel | The approximate minimization of functionals[END_REF] that

λ n -λ L 2 (0,T ;H 1 0 (0,1)) ≤ 2 ν(A r ) ν(A r ) -1 ν(A r ) + 1 n λ 0 -λ L 2 (0,T ;H 1 0 (0,1)) , ∀n ≥ 1 where λ minimizes J ⋆⋆ . ν(A r ) = A r A -1 r
denotes the condition number of the operator A r . Eventually, once the above algorithm has converged we can compute ϕ ∈ Φ as solution of

a r (ϕ, ϕ) + b(ϕ, λ) = l(ϕ), ∀ϕ ∈ Φ.
4 Numerical approximation and experiments 4.1 Some domains q T Let us first describe the domains q T ⊂ Q T := (0, 1) × (0, T ) in which the control is supported we shall use in our numerical experiments. Let γ i T : (0, T ) → (0, 1) be four C ∞ functions defined as follows :

γ 0 T (t) = 3 10 , t ∈ (0, T ), (33)
γ 1 T (t) = 1 2 + 1 10 cos πt T , t ∈ (0, T ), (34)
γ 2 T (t) = β -α T t + α, t ∈ (0, T ), (35)
γ 3 T (t) = 1 2 + 1 4 cos 8πt T , t ∈ (0, T ). ( 36
)
In what follows we choose in (35) α = 0.2 and β = 0.8.

For each i ∈ {0, 1, 2, 3}, let a i T , b i T : (0, T ) → (0, 1) be two functions defined by

a i T (t) = γ i T (t) -δ 0 , b i T (t) = γ i T (t) + δ 0 , t ∈ (0, T ). ( 37
)
for some δ 0 > 0 small enough. We then define the corresponding domains q i T as follows :

q i T = (x, t) ∈ Q T ; a i T (t) < x < b i T (t), t ∈ (0, T ) , i ∈ {0, 1, 2, 3}. (38)
Remark that, in the definition of a i T and b i T , we may consider time-dependent value for δ 0 . Figure 2 display the domains q i T defined by (38) with the controllability time T = 2.2 and δ 0 = 10 -1 . We easily check, that for any T > 2, these domains satisfy the geometric condition of Proposition 2.1.

Discretization

We now turn to the discretization of the mixed formulation (24) assuming r > 0.

Let then Φ h and M h be two finite dimensional spaces parametrized by the variable h such that

Φ h ⊂ Φ, M h ⊂ L 2 (0, T ; H 1 0 (0, 1)), ∀h > 0.
Then, we can introduce the following approximated problems : find

(ϕ h , λ h ) ∈ Φ h × M h solution of a r (ϕ h , ϕ h ) + b(ϕ h , λ h ) = l(ϕ h ), ∀ϕ h ∈ Φ h b(ϕ h , λ h ) = 0, ∀λ h ∈ M h . ( 39
) 0 1 2 t 0 0.5 1 x γ 0 T (t) q 0 T Q T 0 1 2 t 0 0.5 1 x γ 1 T (t) q 1 T Q T 0 1 2 t 0 0.5 1 x γ 2 T (t) q 2 T Q T 0 1 2 t 0 0.5 1 x γ 3 T (t) q 3 T Q T Figure 2:
The time dependent domains q i T , i ∈ {0, 1, 2, 3} defined by (38).

The well-posedness of this mixed formulation is again a consequence of two properties : the coercivity of the bilinear form a r on the subset

N h (b) = {ϕ h ∈ Φ h ; b(ϕ h , λ h ) = 0 ∀λ h ∈ M h }. Actually, from the relation a r (ϕ, ϕ) ≥ r η ϕ 2 Φ , ∀ϕ ∈ Φ
the form a r is coercive on the full space Φ, and so a fortiori on N h (b) ⊂ Φ h ⊂ Φ. The second property is a discrete inf-sup condition : there exists δ h > 0 such that inf

λ h ∈M h sup ϕ h ∈Φ h b(ϕ h , λ h ) ϕ h Φ h λ h M h ≥ δ h . (40)
For any fixed h, the spaces M h and Φ h are of finite dimension so that the infimum and supremum in (40) are reached: moreover, from the property of the bilinear form a r , it is standard to prove that δ h is strictly positive (see Section 4.5). Consequently, for any fixed h > 0, there exists a unique couple (ϕ h , λ h ) solution of (39). On the other hand, the property inf h δ h > 0 is in general difficult to prove and depends strongly on the choice made for the approximated spaces M h and Φ h . We shall analyze numerically this property in Section 4.5.

The finite dimensional and conformal space Φ h must be chosen such that Lϕ h belongs to L 2 (0, T, H -1 (0, 1)) for any ϕ h ∈ Φ h . This is guaranteed for instance as soon as ϕ h possesses second-order derivatives in L 2 loc (Q T ). Therefore, a conformal approximation based on standard triangulation of Q T requires spaces of functions continuously differentiable with respect to both variables x and t.

We introduce a triangulation T h such that Q T = ∪ K∈T h K and we assume that {T h } h>0 is a regular family. We note h := max{diam(K), K ∈ T h } where diam(K) denotes the diameter of K. Then, we introduce the space Φ h as follows :

Φ h = {ϕ h ∈ Φ h ∈ C 1 (Q T ) : ϕ h | K ∈ P(K) ∀K ∈ T h , ϕ h = 0 on Σ T }
where P(K) denotes an appropriate space of polynomial functions in x and t. In this work, we consider for P(K) the reduced Hsieh-Clough-Tocher (HCT for short) C 1 -element. This is a so-called composite finite element and involves 9 degrees of freedom, namely the values of ϕ h , ϕ h,x , ϕ h,t on the three vertices of each triangle K. We refer to [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] page 356 and to [START_REF] Bernadou | Basis functions for general Hsieh-Clough-Tocher triangles, complete or reduced[END_REF][START_REF] Meyer | A simplified calculation of reduced hct-basis functions in a finite element context[END_REF] where the implementation is discussed. We also define the finite dimensional space

M h = {λ h ∈ C 0 (Q T ), λ h | K ∈ Q(K) ∀K ∈ T h , λ h = 0 on Σ T }
where Q(K) denotes the space of affine functions both in x and t on the element K. For any h > 0, we have Φ h ⊂ Φ and M h ⊂ L 2 (0, T ; H 1 0 (0, 1)). For each combination of domains (q T , Q T ) described in Section 4.1 we consider six levels of triangulations T h (numbered from ♯0 to ♯5, from coarser to finer). The number of triangles for some examples of domains q T which will be used the experiments are summarized in Table 1. In Figure 3 we display the meshes ♯1 corresponding to geometries described in Figure 2. In contrast to [START_REF] Cîndea | A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations[END_REF] where the boundary controllability is considered with the constraint Lϕ = 0 as an L 2 (Q T ) function, the equality Lϕ = 0 in Φ is assumed in the weaker space L 2 (0, T ; H -1 (0, 1)). It is not straightforward to handle numerically the scalar product over H -1 which appears in the mixed formulation (39). However, at the finite dimensional level of the discretization, since all the norms are equivalent, a classical trick (see for instance [START_REF] Bochev | Least-squares finite element methods[END_REF][START_REF] Bramble | Least-squares methods for Stokes equations based on a discrete minus one inner product[END_REF]) consists in replacing, for any fixed h, the norm Lϕ h L 2 (0,T ;H -1 (0,1)) by the norm Lϕ h L 2 (Q T ) , up to a constant.

♯ Mesh 0 1 2 3 4 5 q 0 T =2.
In order to do that, first remark that if there exist two constants C 0 > 0 and α > 0 such that

ψ h 2 L 2 (Q T ) ≥ C 0 h α ψ h 2 L 2 (0,T ;H 1 0 (0,1)) , ∀ψ h ∈ Φ h (41)
then a similar inequality it holds for weaker norms. More precisely, we have

ϕ h 2 L 2 (0,T ;H -1 (0,1)) ≥ C 0 h α ϕ h 2 L 2 (Q T ) , ∀ϕ h ∈ Φ h . (42)
Indeed, to obtain (42) it suffices to take ψ h (•, t) = (-∆)

1 2 ϕ h (•, t) in (41). That gives T 0 (-∆) -1 2 ϕ h (•, t) 2 L 2 (0,1) dt ≥ C 0 h α T 0 (-∆) -1 2 ϕ h,x (•, t) 2 L 2 (0,1)
dt.

Since -∆ is a self-adjoint positive operator and ϕ h ∈ Φ h ⊂ H 1 0 (Q T ) we can integrate by parts in both hand-sides of the above inequality and hence we deduce estimate (42). We highlight that the term C 0 h α (and so C 0 and α) does not depend on T .

Assuming that (41) (and consequently (42)) holds (the constants C 0 , α > 0 will be approximated numerically in Section 4.4), we may consider, for any fixed h > 0, the following equivalent definitions of the form a r,h and b h over the finite dimensional spaces Φ h × Φ h and Φ h × M h respectively :

a r,h : Φ h × Φ h → R, a r,h (ϕ h , ϕ h ) = a(ϕ h , ϕ h ) + rC 0 h α Q T Lϕ h Lϕ h dxdt (43) b h : Φ h × M h → R, b h (ϕ h , λ h ) = Q T Lϕ h λ h dxdt. (44)
Let

n h = dim Φ h , m h = dim M h and let the real matrices A r,h ∈ R n h ,n h , B h ∈ R m h ,n h , J h ∈ R m h ,m h and L h ∈ R n h ,1 be defined by a r,h (ϕ h , ϕ h ) = A r,h {ϕ h }, {ϕ h } R n h ,R n h , ∀ϕ h , ϕ h ∈ Φ h , ( 45
) b h (ϕ h , λ h ) = B h {ϕ h }, {λ h } R m h ,R m h , ∀ϕ h ∈ Φ h , ∀λ h ∈ M h , (46)
Q T λ h λ h dx dt = J h {λ h }, {λ h } R m h ,R m h , ∀λ h , λ h ∈ M h , (47)
l(ϕ h ) = L h , {ϕ h } , ∀ϕ h ∈ Φ h ( 48
)
where {ϕ h } ∈ R n h ,1 denotes the vector associated to ϕ h and •, • R n h ,R n h the usual scalar product over R n h . With these notations, the problem (39) reads as follows : find

{ϕ h } ∈ R n h ,1 and {λ h } ∈ R m h ,1 such that A r,h B T h B h 0 R n h +m h ,n h +m h {ϕ h } {λ h } R n h +m h ,1 = L h 0 R n h +m h ,1 . (49)
The matrix A r,h as well as the mass matrix J h are symmetric and positive definite for any h > 0 and any r > 0. On the other hand, the main matrix of order m h + n h in (49) is symmetric but not positive definite. We use exact integration methods developed in [START_REF] Dunavant | High degree efficient symmetrical Gaussian quadrature rules for the triangle[END_REF] for the evaluation of the coefficients of the matrices. The system (49) is solved using the direct LU decomposition method. Let us also mention that for r = 0, although the formulation ( 24) is well-posed, numerically, the corresponding matrix A 0,h is not invertible. In the sequel, we shall consider strictly positive values for r.

Once the approximation ϕ h is obtained, an approximation v h of the control v is given by v h = -ϕ h 1 q T ∈ L 2 (Q T ). The corresponding controlled state y h may be obtained by solving [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF] with standard forward approximation (we refer to [START_REF] Cîndea | Numerical controllability of the wave equation through primal methods and carleman estimates[END_REF], Section 4 where this is detailed). Here, since the controlled state is directly given by the multiplier λ, we simply use λ h as an approximation of y and we do not report here the computation of y h .

Numerical approximation of C 0 and α in (42).

In order to approximate the values of the constants C 0 , α appearing in (41)-(42) we consider the following problem : find α > 0 and C 0 > 0 such that sup

ϕ h ∈Φ h ϕ h 2 L 2 (0,T ;H 1 0 (0,1)) ϕ h 2 L 2 (Q T ) ≤ 1 C 0 h α , ∀h > 0. ( 50
)
Since Φ h is a finite dimensional space, the supremum is, for any fixed h > 0, the solution of the following eigenvalue problem :

∀h > 0, γ h = sup γ : K h {ψ h } = γJ h {ψ h }, ∀{ψ h } ∈ R m h \ {0} (51)
where K h ∈ R n h ,n h and J h ∈ R n h ,n h are the matrices defined by

K h {ψ h }, {ψ h } R n h ,R n h = Q T ψ h,x ψ h,x dxdt, ∀ψ h , ψ h ∈ Φ h , J h {ψ h }, {ψ h } R n h ,R n h = Q T ψ h ψ h dxdt, ∀ψ h , ψ h ∈ Φ h .
We then can choose C 0 and α in (50) such that C 0 h α = γ -1 h , where γ h solves the problem (51). Figure 4 displays γ -1 h corresponding to the matrices K h and J h associated to the domains Q T and q 0 T for the six levels of mesh and T = 2.2. The values of constants C 0 and α which provide the best fitting are C 0 ≈ 1.48 × 10 -2 and α = 2.1993. As expected, we also check that the constant γ h (and so C 0 and α) does not depend on T nor on the controllability domain. From now on, we use these numerical values in the bilinear form a r,h defined by (43).

The discrete inf-sup test

In order to solve the mixed formulation (39), we first test numerically the discrete inf-sup condition (40). Taking η = r > 0 in [START_REF] Cîndea | A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations[END_REF] so that a r,h (ϕ, ϕ) = (ϕ, ϕ) Φ for all ϕ, ϕ ∈ Φ, it is readily seen (see for instance [START_REF] Chapelle | The inf-sup test[END_REF]) that the discrete inf-sup constant satisfies

δ h := inf √ δ : B h A -1 r,h B T h {λ h } = δ J h {λ h }, ∀ {λ h } ∈ R m h \ {0} . (52)
As in the case of boundary controls (see [START_REF] Cîndea | A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations[END_REF]), the matrix B h A -1 r,h B T h is symmetric and positive definite so that the real δ h defined in term of the (generalized) eigenvalue problem (52) is, for any fixed value of the discretization parameter h, strictly positive. This eigenvalue problem is solved using the power iteration algorithm (assuming that the lowest eigenvalue is simple): for any

{v 0 h } ∈ R n h such that {v 0 h } 2 = 1, compute for any n ≥ 0, {ϕ n h } ∈ R n h , {λ n h } ∈ R m h and {v n+1 h } ∈ R m h iteratively as follows : A r,h {ϕ n h } + B T h {λ n h } = 0 B h {ϕ n h } = -J h {v n h } , {v n+1 h } = {λ n h } {λ n h } 2 .
The scalar δ h defined by ( 52) is then given by : 2 reports the values of δ h for various mesh sizes h, for r = 10 -1 and r = 10 -3 and for q T = q 2 2.2 . As expected, we check that δ h decreases as h → 0 and increases as r → 0. More importantly, this table suggests that the sequence δ h remains uniformly bounded by below with respect to h. This property remains true for other control domains q T , as emphasized by Figure 5. 2: δ h vs. h for q T = q 2 2 , r = 10 -1 and r = 10 3 .

δ h = lim n→∞ ( {λ n h } 2 ) -1/2 . Table
♯
We may conclude that the finite elements we use do "pass" the discrete inf-sup test. As we shall see in the next section, this fact implies the convergence of the sequence ϕ h and λ h . • q 0 T with T = 2.2;

+ q 0 T with T = 2;

◭ q 1 T with T = 2; q 2 T with T = 2.2;

q 2 T with T = 2.
◮ q 3

T with T = 2.2.

Figure 5: Values of δ h vs. h for different control domains q i T and r = 10 -1 .

4.6

Numerical experiments for q T = q 2 2 and comparison with the explicit solution

We first consider the domain q T = q 2 2 (see Figure 2) corresponding to an oblique band of length 2δ 0 = 0.2 and T = 2. We define also the following three initial data in V := H 1 0 (0, 1) × L 2 (0, 1):

(EX1) y 0 (x) = sin(πx), y 1 (x) = 0, x ∈ (0, 1), (EX2) y 0 (x) = e -500(x-0.8) 2 , y 1 (x) = 0, x ∈ (0, 1),

(EX3) y 0 (x) = x θ 1 (0,θ) (x) + 1 -x 1 -θ 1 (θ,1) (x), y 1 (x) = 0, θ ∈ (0, 1), x ∈ (0, 1).
In the case where the domain q T depends on the variable t, there is no in general exact solution of the mixed formulation [START_REF] Münch | Optimal distribution of the internal null control for the onedimensional heat equation[END_REF]. However, we can obtain a semi-explicit representation (using Fourier decomposition) of the minimizer (ϕ 0 , ϕ 1 ) of the conjugate functional J ⋆ (see ( 8)), and consequently of the corresponding adjoint variable ϕ, the control of minimal square integrable norm v = -ϕ 1 q T and finally the controlled state y solution of (1-3). In practice, the obtention of the Fourier representation amounts to solve a symmetric linear system. We refer to the Appendix for the details. This allows to evaluate precisely the error vv h L 2 (q T ) with respect to h and confirm the relevance of the method.

Table 3 and 4 collects some numerical values for r = 10 -1 and r = 10 3 respectively corresponding to the initial data (EX1). In the Tables, κ denotes the condition number associated to the linear system (49), independent of the initial data (y 0 , y 1 ). The convergence of vv h L 2 (q T ) , Lϕ h L 2 (0,T ;H -1 (0,1)) and yλ h L 2 (q T ) toward zero as h ց 0 is clearly observed. This is fully in agreement with the uniform discrete inf-sup property we have observed in Section 4.5. We obtain the following rates of convergence with respect to h for r = 10 -1 and r = 10 3 respectively :

r = 10 -1 : v -v h L 2 (q T ) ≈ O(h 1.3 ), Lϕ h L 2 (0,T ;H -1 (0,1)) ≈ O(h 1.3 ), y -λ h L 2 (Q T ) ≈ O(h 1.94 ) r = 10 3 : v -v h L 2 (q T ) ≈ O(h 1.09 ), Lϕ h L 2 (Q T ) ≈ O(h 1.04 ), y -λ h L 2 (Q T ) ≈ O(h 2.01 ).
We refer to Figure 6 which highlights for r = 10 -1 the polynomial convergence of the sequences yλ h L 2 (Q T ) (" ") and vv h L 2 (q 2 T ) ("•") with respect to h. The previous rates suggests that the value of the parameter r has a restricted influence.

♯ Mesh 1 2 3 4 5 h 7.18 × 10 -2 3.59 × 10 -2 1.79 × 10 -2 8.97 × 10 -3 4.49 × 10 -3 v h L 2 (q T ) 5.370 5.047 4.893 4.815 4.776 Lϕ h L 2 (0,T ;H -1 (0,1)) 2.286 9.43 × 10 -1 3.76 × 10 -1 1.5 × 10 -1 6.15 × 10 -2 v -v h L 2 (q T ) 2.45 × 10 -1 9.65 × 10 -2 4.32 × 10 -2 2.29 × 10 -2 1.10 × 10 -2 y -λ h L 2 (Q T ) 5.63 × 10 -3
; q T = q 2 2 ; r = 10 -1 .
The convergence of the method is also observed for the initial data (EX2), mainly supported around x = 0.8 and the less regular data (EX3). Table 5 collects numerical values associated to (EX2), q T = q 2 2 and r = 10 -1 . We obtain the following rates :

vv h L 2 (q T ) ≈ e 5.85 h 1.4 , Lϕ h L 2 (Q T ) ≈ e 7.96 h 1.31 , yλ h L 2 (Q T ) ≈ e 1.508 h 1.62 Figure 7 displays other Q T the dual variable ϕ h and the primal variable λ h for q T = q 2 2.2 . The figures are obtained with the mesh ♯3. As expected, theses variable are mainly concentrated along the characteristics starting from x = 0.8. Table 4: Example EX1; q T = q 2 2 ; r = 10 3 . ; q T = q 2 2.2 ; Norms vv h L 2 (q T ) (•) and yλ h L 2 (Q T ) ( ) vs. h. ; q T = q 2 2.2 : Functions ϕ h (Left) and λ h (Right) over Q T .

♯

Similarly, Table 6 gives the value corresponding to the third example EX3, here with θ = 1/2. We obtain vv h L 2 (q T ) ≈ e 1.69 h 0.53 , Lϕ h L 2 (Q T ) ≈ e 2.88 h 0.56 , yλ h L 2 (Q T ) ≈ e -1.41 h 1.32 . Table 7 gives the numerical results for the Example EX3 with θ = 1/3. We get vv h L 2 (q T ) ≈ e 1.54 h 0.47 , Lϕ h L 2 (Q T ) ≈ e 2.91 h 0.54 , yλ h L 2 (Q T ) ≈ e -1.52 h 1.29 .

Figure 8 displays the dual variable ϕ h and the primal variable λ h for q T = q 2 2.2 and EX3 with θ = 1/3. The figures are again plotted with the mesh ♯3.

v -v h L 2 (q T )
1.4382 8.73 × 10 -1 6.24 × 10 -1 4.24 × 10 -1 3.25 × 10 -1 yλ h L 2 (Q T )

6.86 × 10 -3 3.55 × 10 -3 1.19 × 10 -3 4.64 × 10 -4 1.96 × 10 -4

Table 6: Example EX3 with θ = 1/2; r = 10 -1 ; q T = q 2 2 . ; q T = q 2 2.2 : Functions ϕ h (Left) and λ h (Right).

4.7

Comparison of v L 2 (q T ) for various domains q T with the same measure

The optimization of the support domain q T is particularly relevant in the time dependent situation. As a first step in this direction, we compare numerically in this section the L 2 (q T )-norm of the control v h for various domain q T having the same measure. Along this section, we take r = 10 -1 and T = 2.2. The four domains we consider are q i T =2.2 for i = 0, 1, 2, 3 and are described in Section 4.1.

Table 8 reports the L 2 norms of v h = -ϕ h 1 q i T obtained with the finer mesh (mesh ♯5, see Table 1) associated to each domain.

Figure 9 displays the dual variable ϕ h and the primal one λ h associated to the initial data EX2 and control domains q 3 T . Figure 10 displays the dual variable ϕ h and the primal one λ h associated to the initial data EX3, θ = 1/3 and control domains q 3 T . We remark that any of these domains provides minimal norm controls for every initial data EX1-EX3. In fact, we suspect that the domains minimizing the L 2 -norm of the control of minimal L 2 -norm are strongly connected with the set generated by the characteristics of the initial data. This questions will be investigated in a future study.

v -v h L 2 (q T )
1.3571 9.78 × 10 -1 6.91 × 10 -1 5.13 × 10 -1 3.69 × 10 -1 yλ h L 2 (Q T )

7.12 × 10 -3 3.23 × 10 -3 1.19 × 10 -3 4.82 × 10 -4 2.12 × 10 -4

Table 7: Example EX3 with θ = 1/3; r = 10 -1 ; q T = q 2 2 .

Initial data ) for q T = q i 2.2 , i ∈ {0, 3} for initial data EX1-EX3. 4.8 Behavior of the control as δ 0 ց 0

q 0 T q 1 T q 2 T q 3 T EX1 4
The approach we have developed is valid for any support q T satisfying the hypothesis of Proposition 2.1, in particular arbitrarily thin domain. In this Section we study numerically the evolution of the norm of the controls of minimal L 2 -norm supported in a time dependent domain q T when the measure of these domains goes to 0. Precisely, we consider smaller and smaller values to the parameter δ 0 defining the "thickness" of the domains q T as specified by (37)-(38).

In Table 9 we give the L 2 and L 2 (H -1 ) norms of the controls obtained for the initial data EX1 and control domains q 2 T =2 for δ 0 = 10 -1 2 i for values of i ∈ {0, 1, 2, . . . , 6}. As δ 0 ց 0, q 2 T degenerates into the C 1 -curve γ 2 T defined by (35); this curve satisfies the condition |(γ 2 T ) ′ (t)| < 1 for all t ∈ [0, T ] assumed in [START_REF] Castro | Exact controllability of the 1-D wave equation from a moving interior point[END_REF] to prove the existence of H -1 (∪ t∈(0,T ) γ(t) × {t}) controls.

The numerical values suggest that both norms of the controls are not uniformly bounded (by above) with respect to δ 0 ; this indicates that the L 2 -controllability of (1) with control supported on the curve γ 2 T =2 (see 35) does not hold. Similar behaviors are obtained for the other domains considered in Section 4.1 when δ 0 ց 0. This does not contradict the result of [START_REF] Castro | Exact controllability of the 1-D wave equation from a moving interior point[END_REF] where the H -1 (∪ t∈(0,T ) γ(t) × {t})-controllability is proved in the limit situation. δ 0 10 -1 10 -1 /2 10 -1 /2 2 10 -1 /2 3 10 -1 /2 4 10 -1 /2 9: Example EX1; q T = q 2 2 ; Norms of the control v h obtained for the EX1 for control domains q 2 2 for different values of δ 0 .

Other cases : Non constant velocity of propagation and T < 2

In order to illustrate our approach in a more challenging case we consider the wave equation with a non-constant velocity of propagation c and control supported in a time dependent domain :

   y tt -(c(x)y x ) x = v 1 q T , (x, t) ∈ Q T y(x, t) = 0, (x, t) ∈ Σ T y(x, 0) = y 0 (x), y t (x, 0) = y 1 (x),
x ∈ (0, 1).

(53)

We take the velocity c ∈ C ∞ (0, 1) given by

c(x) =    1, x ∈ [0, 0.45] ∈ [1, 5], (c ′ (x) > 0),
x ∈ (0.45, 0.55) 5,

x ∈ [0.55, 1].

(54)

Note that the Fourier expansion developed in the appendix does not apply in this case. Although the inequality [START_REF] Cui | Exact controllability for a one-dimensional wave equation in non-cylindrical domains[END_REF] is open in this more general case, we observe that the solution of the mixed formulation ( 24) still provides convergent approximations {v h } of controls. Figure 11 depicts the dual variable ϕ h and the primal variable λ h corresponding to the approximation of the control for problem (53), for initial data given by EX3 with θ = 1/3 and control domain q 2 2 . The augmentation parameter is r = 10 -1 .

Since the control acts in a time dependent domain, the geometric controllability condition can hold for values of the controllability time T which are smaller than 2 (we refer to [START_REF] Khapalov | Controllability of the wave equation with moving point control[END_REF]). Figure 12 displays ϕ h and λ h corresponding to the exmaple EX3 for θ = 1/3, T = 1 and q T = q 2 1 . We mention that in this section the domains Q T are discretized using uniform meshes formed by triangles of size h ≈ 10 -2 .

Another, even more challenging situation is the approximation of controls for problem (53) for shorter controllability times. In Figure 13 we display the results obtained for the initial data EX3, domain q 2 T =1 and the velocity of propagation is non constant in space and given by (54). Analyzing the evolution of the norm of λ h with respect to the time, in all the three examples considered in these section it seems to have the controllability, although the hypotheses of Proposition 2.1 are not completely fulfilled.

Conjugate gradient for J ⋆⋆

We illustrate here the Section 3.2: we minimize the functional J ⋆⋆ : L 2 (Q T ) → R with respect to the variable λ. We recall that this minimization corresponds exactly to the resolution of the mixed formulation ( 24) by an iterative Uzawa type procedure. The conjugate gradient algorithm is given at the end of Section 3.2. In practice, each iteration amounts to solve a linear system involving the matrix A r,h of size n h = 4m h (see (49)) which is sparse, symmetric and positive definite. We use the Cholesky method.

We consider the singular situation given by the example EX3 with θ = 0.3, T = 2 and q T = q 2 2 . We take ε = 10 -10 as a stopping threshold for the algorithm (that is the algorithm is stopped as soon as the norm of the residue g n at the iterate n satisfies g n L 2 (0,T ;H 1 0 (0,1)) ≤ 10 -10 g 0 L 2 (0,T ;H 1 0 (0,1)) ) or as the number of iterations is greater than 1 000. The algorithm is initiated with λ 0 = 0 in Q T . Table 10 and 11 display the results for r = 10 -1 and r = 10 3 .

We first check that this iterative method gives exactly the same approximation λ h than the previous direct method (where (49) is solved directly) since, from Proposition (3.1) problem [START_REF] Münch | Optimal distribution of the internal null control for the onedimensional heat equation[END_REF] coincides with the minimization of J ⋆⋆ for r > 0. Then, we observe that the number of iterates is sub-linear with respect to the dimension m h = card({λ h }) of the approximated problem. Once again, this is in contrast with the behavior of the conjugate gradient algorithm when this latter is used to minimize J ⋆ with respect to (ϕ 0 , ϕ 1 ) (see [START_REF] Münch | A uniformly controllable and implicit scheme for the 1-D wave equation[END_REF]).

Figure 14 displays the evolution of the residue g n L 2 (0,T ;H 1 0 (0,1)) / g 0 L 2 (0,T ;H 1 0 (0,1)) with respect to the iteration n for two values of the augmentation parameter : r = 10 -1 and r = 10 3 . The computation has been done with the level mesh ♯3. As expected, we check that a larger value of r improves significantly the convergence of the algorithm: recall that the gradient of J ⋆⋆ in L 2 (H 1 ) is given by: ∇J ⋆⋆ (λ) = A r λ -∆ -1 (Lϕ 0 ) := -∆ -1 (Lϕ) and that r acts on the term Lϕ L 2 (H -1 ) . For a fixe level of mesh, we observe however a lower error λ hy L 2 (Q T ) for r = 10 EX3 with θ = 1/3, for control domain q 2 2 and r = 10 3 .

Concluding remarks and perspectives

We have extended in this work the contribution [START_REF] Cîndea | A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations[END_REF] to a non-cylindrical situation where the support of the controls depend on the time variable. The numerical approximation is based on a direct resolution of the controllability problem through a mixed formulation involving the dual adjoint variable and a Lagrange multiplier, which turns out to coincide with the primal state of the wave equation to be controlled. The well-posedness of this mixed formulation is the consequence of a generalized observability inequality deduced from [START_REF] Castro | Exact controllability of the 1-D wave equation from a moving interior point[END_REF] (and equivalent to the controllability of the equation). The approach leads to a variational formulation over time-space functional Hilbert space without distinction between the time and the space variable and is very appropriate to non-cylindrical situations. At the practical level, the discrete mixed time-space formulation is solved in a systematic way in the framework of the finite element theory: in contrast to the classical approach, there is no need to take care of the time discretization nor of the stability of the resulting scheme, which is often a delicate issue. The resolution amounts to solve a sparse symmetric linear system. As discussed in [START_REF] Cîndea | A mixed formulation for the direct approximation of the control of minimal L 2 -norm for linear type wave equations[END_REF], Section 4.3 (but not employed here), the space-time discretization of the domain allows an adaptation of the mesh so as to reduce the computational cost and capture the main features of the solution.

The numerical experiments reported in this work suggest a very good behavior of the approach: the strong convergence of the sequences {v h } h>0 , approximation of the controls of minimal square integrable norm, are clearly observed as the discretization parameter h tends to zero (as the consequence of the uniform inf-sup discrete property).

As briefly discussed in Section 4.7, this work opens now the possibility to optimize the control v of minimal L 2 (q T )-norm with respect the support q T (equivalently in our case, with respect to the curves a and b, see ( 6)) in the spirit of [START_REF]Optimal design of the support of the control for the 2-D wave equation: a numerical method[END_REF][START_REF] Münch | Optimal distribution of the internal null control for the onedimensional heat equation[END_REF][START_REF] Periago | Optimal shape and position of the support for the internal exact control of a string[END_REF]: for any (y 0 , y 1 ) ∈ H, T > 0 and L ∈ (0, 1), the problem reads : inf q T ∈C L v q T L 2 (q T ) , C L = {q T : q T ⊂ Q T , |q T | = L|Q T | and such that (11) holds} where v q T denotes the control of minimal L 2 (q T ) norm for (1) distributed over q T . Eventually, we also mention that this approach which consists in solving directly the optimality conditions of a controllability problem may be employed to solve inverse problems where, for instance, the solution of the wave equation has to be recovered from a partial observation, typically localized on a sub-domain q T of the working domain: actually, the optimality conditions associated to a least-square type functional can be expressed as a mixed formulation very closed to [START_REF] Münch | Optimal distribution of the internal null control for the onedimensional heat equation[END_REF]. These last two issues will be analyzed in a future work.

A Appendix: Fourier expansion of the control of minimal L 2 (q T )-norm

We expand in this appendix in term of Fourier series the control of minimal L 2 (q T )-norm v for (1) and the corresponding controlled solution y. We use these expansions in Section 4.6 to evaluate with respect to h, the error yλ h L 2 (Q T ) and vv h L 2 (q T ) where the sequence (ϕ h , λ h ) h>0

Figure 3 :

 3 Figure 3: Meshes ♯1 associated with the domains q i T =2.2 : i = 0, 1, 2, 3 from left to right.

Figure 4 :

 4 Figure 4: Values of γ -1 h vs. h (•). The line represents C 0 h α for C 0 ≈ 1.48 × 10 -2 and α ≈ 2.1993.

Figure 6 :

 6 Figure 6: Example EX1; r = 10 -1; q T = q 2 2.2 ; Norms vv h L 2 (q T ) (•) and yλ h L 2 (Q T ) ( ) vs. h.

Figure 7 :

 7 Figure 7: Example EX2; r = 10 -1 ; q T = q 2 2.2 : Functions ϕ h (Left) and λ h (Right) over Q T .

Figure 8 :

 8 Figure8: Example EX3 with θ = 1/3; r = 10 -1 ; q T = q 2 2.2 : Functions ϕ h (Left) and λ h (Right).

Figure 9 : 2 . 2 -Figure 10 : 2 . 2 -

 9221022 Figure 9: Example EX2: q T = q 3 2.2 -Function ϕ h (Left) and λ h (Right) over Q T .

Figure 11 :Figure 12 : 1 -Figure 13 :

 1112113 Figure 11: Example EX3, θ = 1/3:q T = q 22 for a non-constant velocity of propagation -Function ϕ h (Left) and λ h (Right) over Q T .

18 ×

 18 10 -2 3.59 × 10 -2 1.79 × 10 -2 8.97 × 10 -3 4.49 × 10 -3 y L 2 (Q T ) 1.28 × 10 -2 4.77 × 10 -3 1.5 × 10 -3 6.2 × 10 -4 3.52 × 10 -6

r = 10 - 1 Figure 14 :

 10114 Figure 14: Example EX3. Evolution of the residue g n L 2 (0,T ;H 1 0 (0,1)) / g 0 L 2 (0,T ;H 1 0 (0,1)) w.r.t. the iterate n.

Table 1 :

 1 Number of triangles for different meshes and different control domains q i

	2 T =2 q 0 q 1 T =2.2 q 2 T =2.2 q 2 T =2 q 3 T =2.2	207 198 150 179 177 464 1 856 7424 29 696 118 784 475 136 828 3 312 13 248 52 992 211 968 792 3 168 12 672 50 688 202 752 600 2 400 9 600 38 400 153 600 716 2 864 11 456 45 824 183 296 708 2 832 11 328 45 312 181 248

  10 -2 3.59 × 10 -2 1.79 × 10 -2 8.97 × 10 -3 4.49 × 10 -3

	Mesh	1	2	3	4	5
	h 7.18 × r = 10 -1 18.8171	17.5466	17.0642	16.8880	16.8254
	r = 10 3	0.6981	0.8374	0.9246	0.9964	1.0826
		Table

  1.57 × 10 -3 4.04 × 10 -4 1.03 × 10 -4 2.61 × 10 -5 κ 2.46 × 10 7 2.67 × 10 8 2.96 × 10 9 3.03 × 10 10 3.08 × 10 11

Table 3 :

 3 Example EX1

  10 -2 3.59 × 10 -2 1.79 × 10 -2 8.97 × 10 -3 4.49× 10 -3 v h L 2 (q T ) λ h L 2 (Q T ) 9.23 × 10 -2 4.56 × 10 -2 7.70 × 10 -3 1.71 × 10 -3 4.46 × 10 -4 κ 6.12 × 10 8 1.44 × 10 10 1.51 × 10 11 1.55 × 10 12 1.54 × 10 13

	Mesh	1	2	3	4	5
	h	7.18 × 4.1796	4.6185	4.7589	4.7557	4.7291
	Lϕ h L 2 (0,T ;H -1 (0,1))	0.0391	0.0322	0.0162	0.0078	0.0037
	v -v h L 2 (q T ) y	2.4977	1.1341	0.5617	0.2418	0.1201

  10 -2 3.59 × 10 -2 1.79 × 10 -2 8.97 × 10 -3 4.49× 10 -3 v h L 2 (q T ) 10 -1 5.31 × 10 -2 yλ h L 2 (Q T )5.98 × 10 -2 2.78 × 10 -2 8.97 × 10 -3 2.01 × 10 -3 5.38 × 10 -4

	♯ Mesh	1	2	3	4	5
	h	7.18 × 4.8469	7.6514	10.9905	12.6256	12.9022
	Lϕ L 2 (0,T ;H -1 (0,1)) v -v h L 2 (q T )	3.13 × 10 1 8.4949	2.91 × 10 1 6.6975	1.82 × 10 1 3.2515	6.8984 6.24 ×	1.9257

Table 5 :

 5 Example EX2; q T = q 2 2 ; r = 10 -1 . 10 -2 3.59 × 10 -2 1.79 × 10 -2 8.97 × 10 -3 4.49× 10 -3 v h L 2 (q T )

	♯ Mesh	1	2	3	4	5
	h	7.18 × 4.807	4.756	4.707	4.689	4.685
	Lϕ h L 2 (0,T ;H -1 (0,1))	3.858	2.965	1.881	1.232	8.61 × 10 -1

  10 -2 3.59 × 10 -2 1.79 × 10 -2 8.97 × 10 -3 4.49 × 10 -3 v h L 2 (q T )

	♯ Mesh	1	2	3	4	5
	h	7.18 × 5.350	5.263	5.195	5.172	5.165
	Lϕ h L 2 (0,T ;H -1 (0,1))	4.230	3.339	2.095	1.382	1.022

Table 8 :

 8 L 2 -norm v h L 2 (q T

		.3677	3.8770	4.4808	5.5967
	EX2	11.9994 12.0973 10.6268 11.2624
	EX3, θ = 1/3 3.9946	4.5026	5.0132	5.0369

Table 10 :

 10 Conjugate gradient algorithm. EX3 with θ = 1/3, for control domain q 2 2 and r = 10 -1 . 10 -2 3.59 × 10 -2 1.79 × 10 -2 8.97 × 10 -3 4.49× 10 -3 y L 2 (Q T ) 1.15 × 10 -1 5.2 × 10 -21.65 × 10 -2 6.03 × 10 -3 2.89 × 10 -3

	♯ Mesh	1	2	3	4	5
	h 7.18 × ♯ iterate 87	105	119	140	166
	λ h

Table 11 :

 11 Conjugate gradient algorithm.

where (C 1p , C 2p ) p>0 are the Fourier's coefficients of the initial data (y 0 , y 1 ) ∈ V .

solves the discrete mixed formulation (39). We use the characterization of the couple (y, v) in term of the adjoint solution ϕ (see 5), unique minimizer in H of J ⋆ defined by [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF].

We first note (a p , b p ) (p>0) the Fourier coefficients in l 2 (N) × h -1 (N) of the minimizer (ϕ 0 , ϕ 1 ) ∈ L 2 (0, 1) × H -1 (0, 1) of ( 8), such that (ϕ 0 (x), ϕ 1 (x)) = p>0 (a p , b p ) sin(pπx). The adjoint state takes the form ϕ(x, t) = p>0 a p cos(pπt) + b p pπ sin(pπt) sin(pπx).

We get

a p a q q T cos(pπt) cos(qπt) sin(pπx) sin(qπx) dx dt

and

The optimality equation associated to the functional J ⋆ (see 8) then reads

and can be rewritten in term of the (a p , b p ) p>0 as follows :

(57) where M q T denotes a symmetric positive definite matrix and F y0,y1 a vector obtained from (55) and ( 56) respectively. The resolution of the infinite dimensional system (57) (reduced to a finite dimension one by truncation of the sums) allows an approximation of the minimizer (ϕ 0 , ϕ 1 ) of J ⋆ (see ( 8)), and then of ϕ, solution both of the boundary value problem (5) and of the mixed formulation [START_REF] Münch | Optimal distribution of the internal null control for the onedimensional heat equation[END_REF].

Finally, we use that the corresponding control is given by v = -ϕ 1 q T . We expand the corresponding controlled solution as y(x, t) = p>0 b p (t) sin(pπx) where, for all p > 0, the function b p solves the ODE


Pierre-Frederic Villard

Nicolas Koenig

Cyril Perrenot

Manuela Perez

Piers Boshier

Toward a Realistic Simulation of Organ Dissection

Keywords: Dissection, Cutting simulation, Task analysis

Whilst laparoscopic surgical simulators are becoming increasingly realistic they can not, as yet, fully replicate the experience of live surgery. In particular tissue dissection in one task that is particularly challenging to replicate. Limitation of current attempts to simulate tissue dissection include: poor visual rendering; over simplification of the task and; unrealistic tissue properties. In an effort to generate a more realistic model of tissue dissection in laparoscopic surgery we propose a novel method based on task analysis. Initially we have chosen to model only the basic geometrics of this task rather than a whole laparoscopic procedure. Preliminary work has led to the development of a real time simulator performing organ dissection with a haptic thread at 1000Hz. A virtual cutting tool, manipulated through a haptic device, in combination with 1D and 2D soft-tissue models accurately replicate the process of laparoscopic tissue dissection.

Introduction

Surgical procedures using laparoscopic and robot-assisted laparoscopic techniques use specialized instruments and cameras to perform surgical procedures through small (typically 1-2cm) skin incisions. By limiting the extent of surgical trauma in this way in is intended that patient recovery is improved. Unlike conventional open surgery laparoscopic techniques are often technically more demanding for the surgeon. As such major efforts have been made to create high fidelity computer simulations to aid surgeon training.

A major part of any surgical procedure is the separation of tissue through careful dissection. This process is often time consuming and requiring of considerable operator skill. The challenges of accurate tissue dissection are increased in laparoscopic surgery due to fulcrum affect experience when using long instruments pivoted on a body cavity wall, diminished haptic feedback and loss of eye-hand-target axis This particular surgical task is one that is particularly difficult to model as such current computer based simulators do not meet the requirements for advanced surgical training. We propose here a new framework to simulate the tissue dissection task that has benefited from close collaboration between computer scientists and surgeons. The proposed model includes a 3D geometric framework that handles and can be dissected in the same way as in vivo tissues. This article outlines the design strategy and methodology for creating such a simulation as well as the initial results from its testing.

Methods & Materials

State of the art analysis

Tissue dissection in laparoscopic procedures is present in various commercial simulators like the LAP Mentor by Simbionix R and the MIST-VR by Mentice R . However, there remain some drawbacks concerning the organ dissection that alter the simulations realism, including: the interpenetration between the thin peritoneum (membranous lining of the abdominal cavity and organs within) and fat tissues; the triangle edges that give a too regular geometry (see Figure 1.a); the rendering of peritoneum and fat tissues as a continuous as opposed to a heterogeneous surface (see Figure 1.b); the automatic disappearance of fat tissue after being incised, and; the non physiological properties of fat tissues such that they can be excessively stretched without tearing or bleeding. The most significant limitation of existing simulators is that the dissection task is often too easy and too fast compared to a reality.

Task analysis

Dissection could be done with different kinds of tools [START_REF] Mishra | Laparoscopic dissection techniques[END_REF]: standard scissors that cut but do not perform haemostasis; monopolar electrosurgery for coagulation often associated with a hook used to incise; thermofusion devices such as Ligasure R that cauterizes and cuts tissues, and; ultrasonic devices such as Harmonic Scalpel R . We chose to model the latter in our framework (See Figure 2). The procedural steps within this task that have been assessed by surgical experts are presented in Table 1.

Peritoneal structure

Fat tissues are often represented as a continuous medium with a uniform behavior. For example, in one study [START_REF] Punak | Fatty tissue in a haptic illustration environment[END_REF] fat tissue is modeled by a mass-spring system, the geometry of which is computed from a Marching Cube algorithm base on data from real patients. The realism of this model is poor as previously explained.

Steps of the task: [START_REF] Mishra | Laparoscopic dissection techniques[END_REF].

Ensure optimum quality of video image. Factors to consider include: illumination; focus; colour balance; horizontal axis, and; adequacy of field of view [START_REF] Punak | Fatty tissue in a haptic illustration environment[END_REF].

Grasp tissues/organ from which fat tissues are to be dissected. If needed, separate two peritonea in contact [START_REF] Meier | Real-time deformable models for surgery simulation: a survey[END_REF].

Identify intended plane of dissection [START_REF]Soft tissue deformation and cutting simulation for the multimodal surgery training[END_REF].

Ensure appropriate opposition of tissues such that force is applied perpendicular to the intended plane of dissection [START_REF] Forest | Removing tetrahedral from a manifold mesh[END_REF].

Initiate dissection:

• Approach intended plane of dissection with the Harmonic scalpel • Slightly pull upward tissues to dissect • Ensure haemostasis. It requires the tool diathermy activation. Tissues may slightly be attired by the tool • Tissues are slightly compressed • Use dissecting vibrating scissors to divide tissues in the intended plane of dissection • As dissection continues maintain appropriate tension across plane of dissection [START_REF] Nienhuys | A surgery simulation supporting cuts and finite element deformation[END_REF].

Ensure that the tips of all are within the field of view at all times. Open, sharp tipped and diathermy instrument have the potential to cause damage to surrounding structures and tissues. [START_REF] Botsch | A robust procedure to eliminate degenerate faces from triangle meshes[END_REF].

Use of diathermy with lead to the production of smoke that can limit the surgeons view. Under such circumstances, dissection should be temporarily halted until smoke clears or can be evacuated.

We chose instead to model the peritoneum as a heterogeneous structure composed of 1D and 2D models. The material is represented either by a surface within which there are holes or a set of lines constituting a web-like structure.

The 2D structure is in the form of an irregular grid, every point in which may be randomly moved in any direction. Spaces have been removed in the grid to create holes. The Catmull-Clark algorithm is then applied to soften the angles. In comparison the 1D structure comprises of a set of lines composed of several points that are moved, predominantly in the orientation of the tissue plane but also slightly orthogonal to the plane. Skinning is then applied to these lines to increase their thickness. a.

b.

Soft-tissue deformation

A comparison of different soft tissue models has been conducted by the authors of a previous study [START_REF] Meier | Real-time deformable models for surgery simulation: a survey[END_REF]. They concluded that the application of these models is often a compromise based on specific needs.

In the case of this study tissue structures are modeled by mass-spring elements because it is fast enough to handle haptic rendering but also accurate enough to have a realistic behavior. Spring and damping may also be later tuned to simulate inter patient variability.

Cutting simulation

A 1D-cutting simulation is performed on the 1D elements of a web-like structure by subdividing the edge touch by the surgical tool into two spring elements. In order to follow the task analysis an initial axial traction is applied to the two new elements during 1s. The two element parts then follow gravitational and spring like retraction. This is achieved by setting the spring initial length to a value slightly inferior to the edge dimension.

2D cutting simulation is still an open problem but current models are becoming more realistic [START_REF]Soft tissue deformation and cutting simulation for the multimodal surgery training[END_REF]. There are currently three solutions to simulate this procedure. (1) Removing the elements along the cutting path (Figure 3.a). This is very simple to implement and very fast, however, this algorithm gives unrealistic results. It requires a very fine mesh and leads to unconvincing results. In addition, it creates mesh singularities [START_REF] Forest | Removing tetrahedral from a manifold mesh[END_REF]. (2) Subdividing the mesh elements (Figure 3.b). The method consists of breaking down the components in contact with the cutting tool into several basic elements with the edges corresponding to the newly created cutting paths. Although realistic, this method is computationally heavy, due to the fact that whist the soft tissue surface is preserved, the node number dramatically increases, which increases the computation time. Another disadvantage of this method is that it can build degenerate triangles (surface to almost zero), which can lead to unstable behavior. (3) Initial displacement of vertices and separation method (Figure 3.c). This method, used by Nienhuys and van der Stappen [START_REF] Nienhuys | A surgery simulation supporting cuts and finite element deformation[END_REF], consists in first moving the closest point to the cutting path in order to separate the elements. As with the subdivision method, this technique has a high level of realism since cutting coincides exactly with the path of the cutting tool. It has the advantage of keeping the size of the mesh and thus allows the simulation to proceed at constant performance. However, it suffers from the same limitation of the creation of degenerate elements, necessitating triangle removal. Ultimately this third technique was adopted in the current simulation in combination with the algorithm developed by Botsch and Kobbelt [START_REF] Botsch | A robust procedure to eliminate degenerate faces from triangle meshes[END_REF], to suppress degenerated elements. As previously stated the behavior during and after cutting is simulating the real behavior described by the task analysis with the initial traction and the elastic recoil.

Results

The framework has been implemented in C++ using SOFA framework [START_REF] Allard | SOFAan Open Source Framework for Medical Simulation[END_REF] for the organ deformation and for the topology changes. OpenHaptics was also used for the haptic rendering. The interaction with the virtual environment is performed using a Phantom Omni by Sensable R . Figure 2 shows the link between the device (Figure 2.a) and its avatar (Figure 2.b) as well as the features of the surgical tool. Its 3D geometry has been modeled using Blender R .

The geometry strategy is intentionally simple in order to focus on the task. The haptic interaction has a acceptable frame rate of 1000 Hz. In order to guaranty no discontinuity with the haptic feedback while cutting, the device is blocked. It simulates the 2-3 seconds needed for the cauterization.

Various steps of the task analysis are simulated: the initial traction along the axis; the time needed for cauterization, and; the natural elastic behavior of tissues following incision. The matter compression is however not yet implemented and the 3D rendering is not optimal.

Conclusions & Discussion

We have presented a framework for organ dissection. It is based on a close collaboration with surgical experts in order to focus on the realism. Novel outcomes of this work include: (1) a task analysis for organ dissection; (2) a new framework for organ dissection with haptic rendering and soft-tissue deformation, and; (3) an original fat tissue geometry design. The dissection process is based on existing cutting models from the literature with extra heuristic behaviors added in order to reproduce the physiological tissue properties. Early result consists of a basic 3D geometry that can be dissected with a tool following an original cutting simulation process.

Future work will include: improvement of the quality of the 3D rendering; investigation of the influence of different patient's parameters (quantity, elasticity, adherence, etc.). Finally, the model itself will be validated as a clinical tool.

Figure 1 .

 1 Figure 1. Fat tissue and peritoneum: (a) commercial simulator. (b) human patient

Figure 2 .

 2 Figure 2. Surgical tool simulation: (a) Phantom Omni haptic device; (b) 3D surgical tool avatar linked to the haptic device with rotation A depending of the stylus orientation and scissor closure B linked to the button.

Figure 3 .

 3 Figure 3. Cutting method: (a) Element removing, (b) Element subdivision, and ;(c) vertice displacement and separation

Corresponding Author: Pierre-Frederic Villard LORIA Campus scientifique -BP

-54506 Vandoeuvreles-Nancy email: pierrefrederic.villard@loria.fr