Ghislain Haine

AND Denis Matignon

Anass Serhani

NUMERICAL ANALYSIS OF A STRUCTURE-PRESERVING SPACE-DISCRETIZATION FOR AN ANISOTROPIC AND HETEROGENEOUS BOUNDARY CONTROLLED N -DIMENSIONAL WAVE EQUATION AS A PORT-HAMILTONIAN SYSTEM

Keywords: Port-Hamiltonian systems, N -dimensional wave equation, finite element method, structure-preserving discretization, numerical analysis

The anisotropic and heterogeneous N -dimensional wave equation, controlled and observed at the boundary, is considered as a port-Hamiltonian system. A recent structure-preserving mixed Galerkin method is applied, leading directly to a finite-dimensional port-Hamiltonian system: its numerical analysis is carried out in a general framework. Optimal choices of mixed finite elements are then proved to reach the best trade-off between the convergence rate and the number of degrees of freedom for the state error. Exta compatibility conditions are identified for the Hamiltonian error to be twice that of the state error, and numerical evidence is provided that some combinations of finite element families meet these conditions. Numerical simulations are performed in 2D to illustrate the main theorems among several choices of classical finite element families. Several test cases are provided, including non-convex domain, anisotropic or heterogeneous cases and absorbing boundary conditions.

Introduction

The present work addresses the numerical analysis of a structure-preserving space-discretization of an N -dimensional wave equation with boundary control in the formalism of port-Hamiltonian systems. Since it is intended to merge several points of view on the same subject, the authors have taken care to be pedagogical in each section, hence trying to talk to several scientific communities. This choice of presentation will certainly lead readers to find some parts related to his/her domain(s) of research unnecessary. Roughly speaking, this paper is intended to port-Hamiltonian specialists, numerical analysts, and scientific computing users.

1.1. Port-Hamiltonian systems. In the last two decades, infinite-dimensional port-Hamiltonian systems (pHs) [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF][START_REF] Rashad | Twenty years of distributed port-Hamiltonian systems: a literature review[END_REF] have proved to be a very accurate way to model and control complex multi-physics open systems. This framework enjoys several advantages, such as a relevant physical meaning and a useful underlying geometrical structure (namely Stokes-Dirac structure). It has to be pointed out that, even if known Partial Differential Equations (PDEs) are often only rewritten in the pHs formalism in general, this powerful tool also allows a direct modelling of physical systems (see for instance [START_REF] Duindam | Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach[END_REF][START_REF] Serhani | Anisotropic heterogeneous n-D heat equation with boundary control and observation: I. Modeling as port-Hamiltonian system[END_REF][START_REF] Altmann | Port-Hamiltonian formulations of poroelastic network models[END_REF]) which proves useful to derive PDEs. Furthermore, it is intrinsically modular: interaction systems (such as fluidstructure interactions [START_REF] Cardoso-Ribeiro | Port-Hamiltonian model of two-dimensional shallow water equations in moving containers[END_REF], heat-wave interactions [START_REF] Haine | Long-time behavior of a coupled heat-wave system using a structure-preserving finite element method[END_REF], plasma in a tokamak [START_REF] Vu | A structured control model for the thermomagneto-hydrodynamics of plasmas in tokamaks[END_REF], etc.) can be described through the interconnection of several subsystems with a port-Hamiltonian structure, leading to a more complex pHs [START_REF] Cervera | Interconnection of port-Hamiltonian systems and composition of Dirac structures[END_REF][START_REF] Kurula | Dirac structures and their composition on Hilbert spaces[END_REF]. It finally leads to a power balance, expressing the variation of the Hamiltonian functional (often chosen as the system total energy), especially via boundary controls and boundary observations. 1.2. Structure-preserving discretisation. A recent topic of research is to provide accurate (space-) discretization methods to preserve this powerful formalism. Roughly speaking, mainly two non-exclusive communities work on the issue of structure-preserving discretization. The first one makes use of exterior calculus, while the other makes use of vector calculus. It is known that the two points of view are well-founded, and several strategies to merge their advantages efficiently have already been proposed for several discretization issues (see e.g. [START_REF] Kotyczka | Numerical Methods for Distributed Parameter Port-Hamiltonian Systems[END_REF] and the many references therein).

In the present work, a method for the preservation of the power balance of the Hamiltonian (encoded in an underlying Stokes-Dirac structure) is studied. In the wide literature, several strategies have been proposed: we can cite e.g. [START_REF] Hiemstra | High order geometric methods with exact conservation properties[END_REF][START_REF] Seslija | Explicit simplicial discretization of distributed-parameter port-Hamiltonian systems[END_REF][START_REF] Kotyczka | Weak form of Stokes-Dirac structures and geometric discretization of port-Hamiltonian systems[END_REF] for geometric discretizations, [START_REF] Lee | Discrete conservation properties for shallow water flows using mixed mimetic spectral elements[END_REF][START_REF] Egger | Structure preserving approximation of dissipative evolution problems[END_REF] for Galerkin methods and [START_REF] Trenchant | Finite differences on staggered grids preserving the port-Hamiltonian structure with application to an acoustic duct[END_REF] for finite differences method. However, some of these strategies seem difficult to carry over to N -dimensional systems or to apply to complex geometries, while others require post-processing to construct the finite-dimensional Dirac structure. Another structure-preserving community works on the preservation of the de Rham cohomology and related decompositions (such as the Hodge-Helmholtz decomposition). This topic is older and finds its origins in problems such as electromagnetism (see e.g. [START_REF] Monk | Finite element methods for Maxwell's equations[END_REF] and references therein). It is often written in the exterior calculus formalism, allowing for more abstraction, hence more generality, for the construction of discrete differential operators: see e.g. [START_REF] Hiptmair | Discrete Hodge operators[END_REF][START_REF] Bochev | Principles of Mimetic Discretizations of Differential Operators[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] Gerritsma | Discrete Lie derivative[END_REF] for theoretical aspects, and [START_REF] Farle | A port-Hamiltonian finite-element formulation for the Maxwell equations[END_REF][START_REF] Farle | A Port-Hamiltonian Finite-Element Formulation for the Transmission Line[END_REF][START_REF] Gerritsma | Mimetic Spectral Element Method for Anisotropic Diffusion[END_REF] for some applications to partial differential equations.

According to these definitions of structure-preserving discretization, a numerical method for port-Hamiltonian systems should be able to take into account the aforementioned continuous properties at the discrete level. Indeed, this would lead to a relevant physical meaning for the computed quantities (without post-processing), together with an obvious manner to distribute the computations thanks to the modularity property: in particular, each sub-system could be reduced through a structure-preserving model reduction [START_REF] Gugercin | Structure-preserving tangential interpolation for model reduction of port-Hamiltonian systems[END_REF][START_REF] Egger | On Structure-Preserving Model Reduction for Damped Wave Propagation in Transport Networks[END_REF][START_REF] Haine | Structure-Preserving Discretization of a Coupled Heat-Wave System, as Interconnected Port-Hamiltonian Systems[END_REF] prior to their interconnections. Furthermore, in the field of automatic control, several methodologies for efficient control or stabilisation rely on the pHs form of the approximate finite-dimensional system [START_REF] Toledo | Observer-based boundary control of distributed port-Hamiltonian systems[END_REF]: this encourages research for efficient structure-preserving methods of infinite-dimensional pHs, and even more those related to boundary-controlled-andobserved PDEs.

A special case of the mixed Galerkin method, called the Partitioned Finite Element Method (PFEM) [START_REF] Cardoso-Ribeiro | A partitioned finite element method for power-preserving discretization of open systems of conservation laws[END_REF], seems to be one of the most adapted scheme to build a mimetic finite-dimensional Dirac structure [START_REF] Rashad | Twenty years of distributed port-Hamiltonian systems: a literature review[END_REF].

In [START_REF] Joly | Variational Methods for Time-Dependent Wave Propagation Problems[END_REF], the numerical method proposed for the spatial discretization of closed hyperbolic systems, based on the primal-dual or dual-primal formulations given in [START_REF] Joly | Variational Methods for Time-Dependent Wave Propagation Problems[END_REF]Eqs. (15) and ( 16)], and making use of an abstract mixed Galerkin method, can be seen as the starting point of the PFEM for closed systems. Indeed, the idea of partitioning the system to choose on which equation an integration by parts should be applied was already mentioned: "the principle is to multiply the two equations [...] by test functions and to integrate over Ω, but the key point this time is to apply integration by parts only for one of the two equations.", see [35, p. 207]. The new difficulty for port-Hamiltonian systems lies in the boundary terms, namely the control and the observation. The present work investigates the issue of accurate approximations at the boundary. Indeed, the numerical analysis of boundary controlled wave-like systems discretized via the Mixed Finite Element Method (MFEM) often makes use of known results on elliptic systems. This is an easy way to obtain convergence rates, but it definitely strengthens the conditions on the finite element families that can be used, for instance by introducing an artificial need of so-called inf-sup condition.

In the case of Dirichlet boundary control, this makes the numerical schemes artificially complicated, requiring the discretization of a lifting operator (i.e. solving an elliptic system at each time step) or the addition of Lagrange multipliers. To the best of our knowledge, numerical analysis without these difficulties has only been performed on particular choices of finite elements for the former case, see e.g. [START_REF] Joly | Variational Methods for Time-Dependent Wave Propagation Problems[END_REF]Remark 6] and for instance [START_REF] Bécache | An Analysis of New Mixed Finite Elements for the Approximation of Wave Propagation Problems[END_REF][START_REF] Bécache | A New Family of Mixed Finite Elements for the Linear Elastodynamic Problem[END_REF], where new families of mixed finite elements are constructed on purpose. The present work extends this result to open dynamical systems for all conforming finite elements and both Neumann and Dirichlet boundary controls, without such an inf-sup condition, neither the discretization of a lifting operator, nor the use of Lagrange multipliers. The less restrictive conforming conditions (for closed systems) have already been stated in [START_REF] Joly | Variational Methods for Time-Dependent Wave Propagation Problems[END_REF]Eq. (31)], and claimed to be required in [10, Section 7.9] for efficient convergence (with usual finite elements); in [START_REF] Lepe | Locking-free finite element method for a bending moment formulation of Timoshenko beams[END_REF] a similar result has already been obtained for the Timoshenko beam in 1-D. Furthermore, it is shown that adding compatibility conditions between the finite element families in order to preserve de Rham cohomology results in a better convergence rate for the discrete Hamiltonian towards the continuous one, see Theorem 4.4.

Statement of the main results

. The objective of this section is to provide an informal statement of the main result. In this work, the aim is to analyse the convergence of the PFEM, applied on the following system, associated to the N -dimensional anisotropic and heterogeneous wave equation

(1) ρ( - → x ) ∂ 2 t w(t, - → x ) -div T ( - → x ) ---→ grad(w(t, - → x )) = 0, ∀ - → x ∈ Ω, t ≥ 0, w(0, - → x ) = w 0 ( - → x ), ∂ t w(0, - → x ) = w 1 ( - → x ), ∀ - → x ∈ Ω,
together with the following collocated boundary control u and boundary observation y

(2) u(t, - → x ) = T ( - → x ) ---→ grad(w(t, - → x )) - → n ( - → x ), ∀ - → x ∈ ∂Ω, t ≥ 0, y(t, - → x ) = ∂ t w(t, - → x ), ∀ - → x ∈ ∂Ω, t ≥ 0.
In these equations

• Ω is an open bounded domain of R N , N = 1, 2, 3, . . . , with Lipschitz boundary ∂Ω; • -→ n is the outward normal at the boundary ∂Ω;

• w(t, -→ x ) is the deflection from the equilibrium position at point -→ x ∈ Ω and time t ≥ 0; • u is the boundary control corresponding to forces applied at the boundary; • y is the collocated boundary observation corresponding to the measured velocities on ∂Ω; • ρ is the mass density, supposed to be bounded from above and below (almost everywhere) by ρ + and ρ -> 0 respectively; • T is Young's elasticity modulus, supposed to be a real symmetric tensor bounded from above and below (almost everywhere, in the matrix-norm sense) by T + I and T -I respectively, where T -> 0 and I is the identity tensor;

• stands for the transpose of vectors or matrices.

We associate to system (1)-( 2) the Hamiltonian

H(t) := 1 2 Ω ---→ grad w(t, - → x ) T ( - → x ) ---→ grad w(t, - → x ) + ρ( - → x ) ∂ t w(t, - → x ) 2 d - → x ,
made of the potential and kinetic energies of the physical system.

The following statement is an abridged and informal version of Theorems 4.2 and 4.4.

Main results. Let us denote the strain -→ α q := ---→ grad(w), the linear momentum α p := ρ∂ t w, and their discrete counterparts -→ α d q and α d p obtained by the Partitioned Finite Element Method. The Hamiltonian then rewrites (3)

H(t) := H( - → α q (t), α p (t)) := 1 2 Ω - → α q (t, - → x ) T ( - → x ) - → α q (t, - → x ) + α p (t, - → x ) 2 ρ( - → x ) d - → x .
Let us define the discrete Hamiltonian

H d (t) := H( - → α d q (t), α d p (t)
) and the two errors:

• E X (t) := - → α q (t) α p (t) - - → α d q (t) α d p (t) X
the absolute error in a suitable energy space X ; • E H (t) := H(t) -H d (t) the error between the continuous and the discrete Hamiltonians.

Under suitable assumptions (regularity, conformity and order given by a parameter κ) on the three finite element families (for -→ α q in Ω, α p in Ω, and (u, y) on ∂Ω), for all T > 0, all initial data smooth enough and all u smooth enough, there exist C > 0, independent of h the mesh size parameter, and h * > 0 such that

E X (t) ≤ C h κ , ∀t ∈ [0, T ], h ∈ (0, h * ).
Furthermore, under compatibility assumptions between the three finite element families, the Hamiltonian error satisfies

(4) E H (t) -E H (0) = 1 2 E X (t) 2 -E X (0) 2 , ∀t ∈ [0, T ],
giving a convergence rate of order 2κ.

1.4. Organization of the paper. The paper is organized as follows: in Section 2, the well-posedness of the physical system (1)-( 2) is recalled, and some assumptions on the regularity of the solutions are made. In Section 3, the PFEM is applied and discussed, and the resulting finite-dimensional Dirac structure is highlighted. The case of Dirichlet boundary control, i.e. switching control and observation, is addressed in Section 3.2. In Section 4, the main convergence results are proved for a general Galerkin approximation method, namely Theorem 4.2 for the state error, and Theorem 4.4 for the Hamiltonian error. In Section 5, accurate combinations of finite elements are proposed for a given rate of convergence, by minimizing the number of degrees of freedom. In Section 6, 2D simulations are provided to exhibit the proven convergence rates and its optimality (i.e. maximizing the convergence rate with the minimal number of degrees of freedom); several test cases are provided to illustrate the flexibility of the method. Finally, Section 7 concludes this work with a summary of the results and draws some perspectives.

The N -dimensional wave equation as a pHs

In this section, the boundary-controlled-and-observed wave system (1)-( 2) is firstly recast as a port-Hamiltonian system. This system has already been studied in the pHs framework in [START_REF] Kurula | Linear wave systems on n-D spatial domains[END_REF], in a more general context, i.e. with several boundary conditions on a partition of ∂Ω and internal fluid damping. Secondly, well-posedness is recalled [START_REF] Kurula | Linear wave systems on n-D spatial domains[END_REF] and a refined regularity result is conjectured, assuming a higher regularity of the physical parameters, the initial data, and the control.

Although it should be possible to prove the regularity assumptions making use of Boundary Control Systems theory [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]Chapter 10.] or by adapting the results in [START_REF] Kurula | Linear wave systems on n-D spatial domains[END_REF] to the uniform boundary control case considered here, it goes beyond the scope of this work.

2.1. The distributed-parameters port-Hamiltonian system. From now on, H κ (Ω) denotes the usual Sobolev space for κ ∈ R, and H 0 (Ω) is identified with L 2 (Ω) (the same notations are used on the boundary ∂Ω). We also write L 2 (Ω) := (L 2 (Ω)) N and H κ (Ω) := (H κ (Ω)) N . Definition 2.1 (Traces [16,Chapter 2]). The linear trace operators are defined as follows

• the Dirichlet trace operator γ 0 , defined by γ 0 (v) := v |∂Ω for v ∈ C ∞ (Ω), extends continuously from H 1 (Ω) onto H 1 2 (∂Ω); • the normal trace operator γ ⊥ , defined by γ ⊥ ( - → v ) := ( - → v - → n ) |∂Ω on (C ∞ (Ω)) N , extends continuously from H(div; Ω) := - → v ∈ L 2 (Ω) | div - → v ∈ L 2 (Ω) , onto H -1 2 (∂Ω).
The so-called Green's formula then reads: for all -

→ v ∈ H(div; Ω), v ∈ H 1 (Ω), (5)
Ω v( - → x )div( - → v ( - → x )) d - → x = - Ω ---→ grad(v( - → x )) - → v ( - → x ) d - → x + γ ⊥ ( - → v ), γ 0 (v) H -1 2 (∂Ω),H 1 2 (∂Ω)
.

The last term in [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF] is the duality bracket between H 1 2 (∂Ω), and H -1 2 (∂Ω). Note that as soon as γ ⊥ ( -→ v ) ∈ L 2 (∂Ω), this bracket reduces to the usual L 2 (∂Ω)-inner product [55, 2.9]. Let us define the strain -→ α q := ---→ grad(w) and the linear momentum α p := ρ ∂ t w. Then one can rewrite the first line of System (1) as

(6) ∂ t - → α q ∂ t α p = 0 ---→ grad div 0 =:J T 0 0 ρ -1 =:Q - → α q α p .
The boundary control and observation (2) then read

(7) u = γ ⊥ T - → α q , y = γ 0 ρ -1 α p .
We already define in (3) the Hamiltonian of ( 6)-( 7)

H(t) := H( - → α q (t), α p (t)) := 1 2 Ω - → α q (t, - → x ) T ( - → x ) - → α q (t, - → x ) + α p (t, - → x ) 2 ρ( - → x ) d - → x ,
corresponding to the sum of the potential and kinetic energy, i.e. the total mechanical energy of the system. Making use of Green's formula [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF] together with ( 6)-( 7), one gets that for all t ≥ 0 (8)

d d t H(t) = u(t), y(t) H -1 2 (∂Ω),H 1 2 (∂Ω)
, meaning that the variation of energy is the power supplied to the system at the boundary [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF].

Remark 2.1. In the Hamiltonian formalism, -→ α are called the energy variables while -→ e := δ-→ α H( -→ α ), the variational derivative of H with respect to -→ α [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF] are the co-energy variables. The relations between -→ e and -→ α , linear in the present case, are known as the constitutive relations, which enable to close the system of equations.

Remark 2.2. The co-energy variables -→ e q := δ-→ α q H and e p := δ αp H are physically meaningful: -→ e q = T -→ α q = T ---→ grad (w) is the stress, and e p = ρ -1 α p = ∂ t w is the deflection velocity. Furthermore, as seen in ( 8), a relevant way to control and observe the system through the boundary is related to both the traces of these co-energy variables.

2.2. Existence and uniqueness of solutions. Let X := L 2 (Ω) × L 2 (Ω) be the energy space, endowed with the inner product

- → z 1 , - → z 2 X := Q - → α q 1 α p 1 , - → α q 2 α p 2 L 2 (Ω)×L 2 (Ω) = Ω T ( - → x ) - → α q 1 ( - → x ) - → α q 2 ( - → x ) + ρ( - → x ) -1 α p 1 ( - → x ) α p 2 ( - → x ) d - → x , for all ( - → z 1 , - → z 2 ) := - → α q 1 α p 1 , - → α q 2 α p 2 ∈ X 2 .
It is clear from the assumption on ρ and T that the norm inherited from this inner product is equivalent to the usual

L 2 (Ω) × L 2 (Ω)-norm.
In addition, let

Z := Q -1 H(div; Ω) H 1 (Ω)
be the solution space, U := H -1 2 (∂Ω) the control space and Y := U = H 1 2 (∂Ω) the observation space. It has been shown in [START_REF] Kurula | Linear wave systems on n-D spatial domains[END_REF]Corollary 4.3] that this leads to an internally well-posed strong impedance conservative boundary control system on (U, X , Y). Theorem 2.1 (Corollary 4.3 in [START_REF] Kurula | Linear wave systems on n-D spatial domains[END_REF]).

For all u ∈ C 2 ([0, ∞); U), - → z 0 := - → α q 0 α p 0 := ---→ grad(w 0 ) ρ -1 w 1 ∈ Z such that u(0) = γ ⊥ T ---→ grad(w 0 )
, there exists a unique solution to (6)- [START_REF] Bécache | An Analysis of New Mixed Finite Elements for the Approximation of Wave Propagation Problems[END_REF] with With this material at hand, we are now able to define properly what is meant by formally skew-symmetric: the operator J Q restricted to H 0 (div; Ω) × H 1 (Ω), where H 0 (div; Ω

- → z = - → α q α p ∈ C 1 ([0, ∞); X ) ∩ C([0, ∞); Z), y ∈ C([0, ∞); Y).
) := - → v ∈ H(div; Ω) | γ ⊥ - → v = 0
, is skew-adjoint on X . This follows from Green's formula [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF].

As is often the case in numerical analysis, sufficient regularity of the solution of the continuous problem is required to be allowed to use the interpolation error inequalities and prove convergence. The following development is formal, and its usefulness for the rest of this work will be enlightened in Remark 2.4.

For all integer κ ≥ 0, assume that ∂Ω is C κ+2 and define

• H 0 (div; Ω) := L 2 (Ω) and

H κ (div; Ω) := - → v ∈ H κ-1 (Ω) | div - → v ∈ H κ-1 (Ω) , when κ ≥ 1, endowed with the inner product - → v 1 , - → v 2 H κ (div;Ω) := - → v 1 , - → v 2 H κ-1 (Ω) + div( - → v 1 ), div( - → v 2 ) H κ-1 (Ω) ; • X κ := H κ (div; Ω) × H κ (Ω) endowed with the following bilinear form, for all - → z i := - → v i v i ∈ X κ , i = 1, 2 - → z 1 , - → z 2 Xκ := Q - → z 1 , - → z 2 H κ (div;Ω)×H κ (Ω) := T - → v 1 , - → v 2 H κ (div;Ω) + ρ -1 v 1 , v 2 H κ (Ω) ,
the energy space;

• Z κ := Q -1 H κ+1 (div; Ω) H κ+1 (Ω) the solution space;

• U κ := H κ-1 2 (∂Ω) the control space; • Y κ := H κ+ 1 2 (∂Ω) the observation space. It is known from [16, Chapter 2; Theorem 1 & Proposition 10] that the traces of Definition 2.1 satisfy

• γ 0 is continuous from H κ+1 (Ω) onto H κ+ 1 2 (∂Ω); • γ ⊥ is continuous from H κ+1 (div; Ω) onto H κ-1 2 (∂Ω).
Assume furthermore that ρ and T are smooth enough for X κ to be a Hilbert space.

Conjecture 2.1. With the above notations and assumptions, it holds

∀ - → z 0 := ---→ grad(w 0 ) ρ -1 w 1 ∈ Z κ , ∀ u ∈ C 2 ([0, ∞); U κ ) : u(0) = γ ⊥ T ---→ grad(w 0 ) ,
there exists a unique solution -→ z satisfying

(H0) - → z ∈ C 1 ([0, ∞); X κ ) ∩ C([0, ∞); Z κ ), and y ∈ C([0, ∞); Y κ ).
Roughly speaking, it claims that increasing the regularity on ρ, T , ∂Ω, w 0 , w 1 , and u, increases the space regularity of solutions (using the continuity and surjectivity of γ 0 and γ ⊥ ). This seems legitimate according to [START_REF] Kurula | Linear wave systems on n-D spatial domains[END_REF]Corollary 4.3] (which includes the above case κ = 0). Although this is not proved, this seems to be a reasonable conjecture at the mathematical level, i.e. for existence and uniqueness of smooth solutions.

Remark 2.4. The main purpose of (H0) is to provide a relation between the maximal H (div; Ω)-and H k (Ω)-regularities of -→ α q and α p respectively, allowing for an optimal choice of the order of the finite-dimensional spaces of approximation.

Finally, it has to be kept in mind that Conjecture 2.1 implies that if κ is supposed to be the maximal regularity (in space) of -→ α q , i.e.

- → α q (t) ∈ T -1 H κ (div; Ω) but - → α q (t) ∈ T -1
H κ+1 (div; Ω), then κ is also the maximal regularity (in space) of α p , i.e. α p (t) ∈ ρH κ (Ω) but α p (t) ∈ ρH κ+1 (Ω), and reciprocally.

The weak co-energy formulation.

As it is intended to apply a conforming finite element method, it is mandatory to use a formulation which enables the use of available finite elements in numerical softwares without destroying the sparsity property of the FEM, e.g. avoiding matrix inversion.

A simple way to achieve this is to work on the co-energy formulation, which transfers the physical parameters from the right-hand side of (6) to its left-hand side, by inverting them at the continuous level. More precisely, rewriting the initial system (6) by making use of the relations -→ α q = T -1 -→ e q and α p = ρe p , one obtains the equivalent port-Hamiltonian system, known as the co-energy formulation

(9) T -1 0 0 ρ ∂ t - → e q ∂ t e p = J - → e q e p , u = γ ⊥ - → e q , y = γ 0 (e p ) .
At the discrete level, this gives rise to weighted mass matrices, carrying all the physical parameters.

Multiplying (9) in L 2 (Ω) × L 2 (Ω) by arbitrary test functions

- → v q v p , one gets T -1 0 0 ρ ∂ t - → e q ∂ t e p , - → v q v p L 2 (Ω)×L 2 (Ω) = 0 ---→ grad div 0 - → e q e p , - → v q v p L 2 (Ω)×L 2 (Ω)
.

which also reads      ∂ t - → e q , T -1 - → v q L 2 (Ω) = ---→ grad (e p ) , - → v q L 2 (Ω) , ∂ t e p , ρ v p L 2 (Ω) = div - → e q , v p L 2 (Ω) .
At this stage, the boundary control in [START_REF] Bochev | Principles of Mimetic Discretizations of Differential Operators[END_REF] does not appear in the formulation yet.

To this end, we apply Green's formula (5) on the second line only and obtain (10)

     ∂ t - → e q , T -1 - → v q L 2 (Ω) = ---→ grad (e p ) , - → v q L 2 (Ω) , ∂ t e p , ρ v p L 2 (Ω) = -- → e q , ---→ grad (v p ) L 2 (Ω) + u, γ 0 (v p ) U ,Y , remembering that U = H -1 2 (∂Ω) and Y = H 1 2 (
∂Ω). These equations make sense if -→ v q ∈ L 2 (Ω) and v p ∈ H 1 (Ω). Altogether, the test functions have to belong to

L 2 (Ω) × H 1 (Ω) (note that this is neither X nor Z defined in Section 2.2).

Structure-preserving discretisation

The aim of this section is to show how a simple integration by part on a partition of the state space, as in the MFEM, is able to transforms a distributed-parameters port-Hamiltonian system (i.e. infinite-dimensional) into a lumped-parameters port-Hamiltonian system (i.e. finite-dimensional). The main difference with the MFEM is the non-homogeneous boundary condition (i.e. the boundary control) applied to the partition to be integrated by part. The major interest of this scheme, known as the Partitioned Finite Element Method (PFEM) [START_REF] Cardoso-Ribeiro | A partitioned finite element method for power-preserving discretization of open systems of conservation laws[END_REF], is that it directly leads to a discrete version of the power-balance satisfied by the discrete Hamiltonian, defined as the continuous one evaluated on the approximations of energy variables.

3.1. Discrete weak formulation and matrix form. We are now in a position to discretize the system in space. Assume that we have at our disposal three finite dimensional spaces, typically given by finite elements respectively L 2 (Ω)conforming, H 1 (Ω)-conforming and H 1 2 (∂Ω)-conforming

H q := Span - → ϕ i q i=1,...,Nq ⊂ L 2 (Ω) of dimension N q ∈ N, H p := Span ϕ k p k=1,...,Np ⊂ H 1 (Ω) of dimension N p ∈ N, and
H ∂ := Span (ψ m ) m=1,...,N ∂ ⊂ H 1 2 (∂Ω) of dimension N ∂ ∈ N.
Remark 3.1. Note that the boundary basis (ψ m ) m=1,...,N ∂ , used to approximate both u and y for simplicity, is chosen to be H 1 2 (∂Ω)-conforming, while u only required H -1 2 (∂Ω) to be approximated in a conforming manner. This is a convenient assumption which leads to the usual L 2 (∂Ω)-inner product at the boundary, without loss of generality as soon as κ ≥ 1.

Let us approximate the function -→ e q in H q by

- → e q (t, - → x ) - → e d q (t, - → x ) := Nq i=1 e i q (t) - → ϕ i q ( - → x ) = - → Φ Φ Φ q ( - → x ) e q (t), ∀t ≥ 0, - → x ∈ Ω,
where for all t ≥ 0 and -→ x ∈ Ω, we introduce the compact notations

- → Φ Φ Φ q ( - → x ) :=      - → ϕ 1 q ( - → x ) . . . - → ϕ Nq q ( - → x )      ∈ R Nq×N , e q (t) :=    e 1 q (t) . . . e Nq q (t)    ∈ R Nq .
In the same way, e p is approximated in H p by

e p (t, - → x ) e d p (t, - → x ) := Np k=1 e k p (t)ϕ k p ( - → x ) = Φ p ( - → x ) e p (t), ∀t ≥ 0, - → x ∈ Ω,
where for all t ≥ 0 and all -→ x ∈ Ω, we have the compact notations

Φ p ( - → x ) :=    ϕ 1 p ( - → x ) . . . ϕ Np p ( - → x )    ∈ R Np , e p (t) :=    e 1 p (t) . . . e Np p (t)    ∈ R Np . Finally, u is approximated in H ∂ by u(t, - → s ) u d (t, - → s ) := N ∂ m=1 u m (t)ψ m ( - → s ) = Ψ( - → s ) u(t), ∀t ≥ 0, - → s ∈ ∂Ω,
where for all t ≥ 0 and all -→ s ∈ ∂Ω, we also use the compact notations

Ψ( - → s ) :=    ψ 1 ( - → s ) . . . ψ N ∂ ( - → s )    ∈ R N ∂ , u(t) :=    u 1 (t) . . . u N ∂ (t)    ∈ R N ∂ .
It is now possible to formulate the discrete variational formulation from the continuous one [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF] on H q × H p × H ∂ : for all j = 1, . . . , N q and all = 1, . . . , N p , we are seeking for ( -

→ e d q , e d p ) ∈ H q × H p such that (11)      ∂ t - → e d q , T -1 - → ϕ j q L 2 (Ω) = ---→ grad e d p , - → ϕ j q L 2 (Ω) , ∂ t e d p , ρ ϕ p L 2 (Ω) = -- → e d q , ---→ grad ϕ p L 2 (Ω) + u d , γ 0 ϕ p L 2 (∂Ω) .
From the definition of -→ e d q , e d p and u d , this leads to

(12)                        Nq i=1 d d t e i q - → ϕ i q , T -1 - → ϕ j q L 2 (Ω) = Np k=1 e k p ---→ grad ϕ k p , - → ϕ j q L 2 (Ω) , Np k=1 d d t e k p ϕ k p , ρ ϕ p L 2 (Ω) = - Nq i=1 e i q - → ϕ i q , ---→ grad ϕ p L 2 (Ω) + N ∂ m=1 u m ψ m , γ 0 ϕ p L 2 (∂Ω) . Now, denoting ---→ grad (Φ p ) :=       ---→ grad ϕ 1 p . . . ---→ grad ϕ Np p       ∈ R Np×N ,
the gradient of the p-type family (which is H 1 (Ω)-conforming by hypothesis),

M T -1 := Ω - → Φ Φ Φ q ( - → x ) T -1 ( - → x ) - → Φ Φ Φ q ( - → x ) d - → x ∈ R Nq×Nq , M ρ := Ω ρ( - → x ) Φ p ( - → x ) Φ p ( - → x ) d - → x ∈ R Np×Np ,
the mass matrices taking the metric of X into account,

D := Ω - → Φ Φ Φ q ( - → x ) ---→ grad Φ p ( - → x ) d - → x ∈ R Nq×Np ,
the averaged gradient and

B ∂ := ∂Ω γ 0 (Φ p ) ( - → s ) Ψ( - → s ) d - → s ∈ R Np×N ∂ ,
the discrete boundary control operator, we get the following finite-dimensional dynamical system from (12) (13)

M T -1 0 0 M ρ d d t e q (t) e p (t) = 0 D -D 0 e q (t) e p (t) + 0 B ∂ u(t), ∀t ≥ 0.
Finally, defining

M ∂ := ∂Ω Ψ( - → s ) (Ψ( - → s )) d - → s ∈ R N ∂ ×N ∂ ,
the boundary mass matrix and

B ∂ := 0 B ∂ ∈ R (Nq+Np)×N ∂ ,
the extended boundary control operator, the output is given for all t ≥ 0 by ( 14)

M ∂ y(t) := B ∂ e q (t) e p (t) = B ∂ e p (t).
Now, system ( 13)-( 14) is a finite-dimensional port-Hamiltonian system.

Remark 3.2. One can gather equations ( 13)-( 14) under the flows-efforts formulation often used in the port-Hamiltonian systems community ( 15)

  M T -1 0 0 0 M ρ 0 0 0 M ∂        d d t e q (t) d d t e p (t) -y(t)      =   0 D 0 -D 0 B ∂ 0 -B ∂ 0     e q (t) e p (t) u(t)   .
The block diagonal symmetric positive-definite matrix constituted by the mass matrices on the left-hand side carry the physical parameters by taking the induced metric into account. On the right-hand side, the skew-symmetric matrix is known as the extended structure matrix. Together, they give a kernel representation [START_REF] Van Der Schaft | Port-Hamiltonian Systems Theory: An Introductory Overview[END_REF] of the underlying Dirac structure [START_REF] Serhani | A Partitioned Finite Element Method for the Structure-Preserving Discretization of Damped Infinite-Dimensional Port-Hamiltonian Systems with Boundary Control[END_REF].

At this stage, the co-energy formulation (9) of the infinite-dimensional port-Hamiltonian system (6)-( 7) has been accurately discretized as a finite-dimensional port-Hamiltonian system (13)- [START_REF] Cardoso-Ribeiro | Port-Hamiltonian model of two-dimensional shallow water equations in moving containers[END_REF]. The next step is to define a discrete version of the Hamiltonian H in order to perfectly mimic the power balance (8). Definition 3.1. The discrete Hamiltonian H d is defined as the evaluation of H on the approximations -→ α d q := T -1 -→ e d q , and α d p := ρ e d p , namely

H d (t) := H( - → α d q (t), α d p (t)). Proposition 3.1. The discrete Hamiltonian H d reads (16) H d (t) = 1 2 e q (t) M T -1 e q (t) + 1 2 e p (t) M ρ e p (t).
For all t ≥ 0, the following power balance holds

(17) d d t H d (t) = (u(t)) M ∂ y(t), = u d (t), y d (t) L 2 (∂Ω) ,
which is the discrete counterpart of (8).

Proof. The equality (16)

u(t, - → x ) = ∂ t w(t, - → x ), ∀ - → x ∈ ∂Ω, t ≥ 0, y(t, - → x ) = T ( - → x ) ---→ grad(w(t, - → x )) - → n ( - → x ), ∀ - → x ∈ ∂Ω, t ≥ 0, leading to (10S)      ∂ t - → e q , T -1 - → v q L 2 (Ω) = -e p , div - → v q L 2 (Ω) + γ ⊥ - → v q , u L 2 (∂Ω) , ∂ t e p , ρ v p L 2 (Ω) = div - → e q , v p L 2 (Ω) ,
instead of [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF]. The PFEM would then provide the following matrices

D = - Ω div - → Φ Φ Φ q ( - → x ) Φ p ( - → x ) d - → x , B ∂ = ∂Ω Ψ( - → s ) γ ⊥ - → Φ Φ Φ q ( - → s ) d - → s , such that   M T -1 0 0 0 M ρ 0 0 0 M ∂        d d t e q (t) d d t e p (t) -y(t)      =    0 D B ∂ -D 0 0 -B ∂ 0 0      e q (t) e p (t) u(t)  
In this case, H q is chosen H(div; Ω)-conforming and H p only L 2 (Ω)-conforming.

Remark 3.3. More complex causalities, such as mixed boundary controls, boundary damping, etc., can be handled in the same manner, as it will be done in Section 6.5 for absorbing boundary condition. We refer to [START_REF] Serhani | Partitioned Finite Element Method for port-Hamiltonian systems with Boundary Damping: Anisotropic Heterogeneous 2-D wave equations[END_REF][START_REF] Brugnoli | Partitioned Finite Element Method for Power-Preserving Structured Discretization with Mixed Boundary Conditions[END_REF] for more details.

Numerical analysis

This section is the core of this work: we state the main theorems and provide their proof. They are given under usual assumptions for Galerkin methods, gathered below from (H1) to (H5), that prove classical for the finite element method (see e.g. [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF][START_REF] Gatica | A Simple Introduction to the Mixed Finite Element Method: Theory and Applications[END_REF]). The aim of such an abstract numerical analysis is to propose a general framework to deal with several kinds of approximation families at the same time. This allows for recent developments such that, e.g. , conforming discontinuous Galerkin elements on rectangular mesh [START_REF] Feng | A conforming discontinuous Galerkin finite element method on rectangular partitions[END_REF], or even for meshfree methods [START_REF] Oliveira | Weak-form collocation -A local meshless method in linear elasticity[END_REF]. Next, Section 5 shall focus on well-known examples of suitable choices of such families of conforming finite elements.

Let us consider the following

• E X (t) := - → α q (t) α p (t) - - → α d q (t) α d p (t) X
the absolute error in X between the continuous and the discrete energy variables;

• E H (t) := H(t) -H d (t) the error between the continuous and the discrete Hamiltonians. The aim is to analyse the asymptotic behaviour of those errors, when the values of N q , N p and N ∂ tend towards ∞ (e.g. when the mesh size parameter tends towards 0). Furthermore, the best trade-off between the discretization orders of H q , H p and H ∂ are provided: the number of degrees of freedom is minimized for each fixed desired convergence rate.

Remark 4.1. Thanks to [START_REF] Duindam | Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach[END_REF], it holds E H (t) := H(t) -H d (t) = H(0) -H d (0) for all t ≥ 0, as soon as the system is closed (i.e. with u ≡ 0). This result is well-known since several decades using the MFEM [START_REF] Joly | Variational Methods for Time-Dependent Wave Propagation Problems[END_REF].

Remark 4.2. Theorem 4.2 could be obtained by considering the homogeneous isotropic case only (corresponding to an identification between energy and co-energy variables). However, anisotropy and heterogeneity induce mandatory modifications about the compatibility conditions allowing for preserving de Rham cohomology, as will be shown in Theorem 4.4.

4.1. Notations, hypotheses and basic properties. In the sequel, the following general hypotheses are assumed. These assumptions are made of

• usual Galerkin estimates on H 1 and L 2 ;

• an inverse inequality between the H 1 -and L 2 -norms on the finite-dimensional space H p ⊂ H 1 (Ω); • an estimate of the L 2 -projection in the H 1 -norm, which proves useful to get optimality.

Notations. Let us denote

• h ∈ (0, h * ) a parameter vow to tends to 0, where h * > 0 (i.e. h is small enough);

• P p the L 2 (Ω)-orthogonal projector from L 2 (Ω) onto H p ; • P 1,p the H 1 (Ω)-orthogonal projector from H 1 (Ω) onto H p ; • - → P q the L 2 (Ω)-orthogonal projector from L 2 (Ω) onto H q .
In order to take into account the metric induced by the operator Q on X , we also introduce • -→ P P P q the orthogonal projector from L 2 (Ω) endowed with the weighted inner

product - → v 1 , T - → v 2 L 2 for all - → v 1 , - → v 2 ∈ L 2 (Ω), onto V q := T -1
H q ;

• P p the orthogonal projector from L 2 (Ω) endowed with the weighted inner product

v 1 , ρ -1 v 2 L 2 for all v 1 , v 2 ∈ L 2 (Ω), onto V p := ρH p .

Hypotheses.

There exists h * > 0 such that for all κ ≥ 0, (H1) ∃C p > 0, ∃θ p ≥ 0 :

P p v p -v p L 2 (Ω) ≤ C p h θp v p H κ+1 (Ω) , ∀v p ∈ H κ+1 (Ω), ∀h ∈ (0, h * ); (H2) ∃C 1,p > 0, ∃θ 1,p ≥ 0 : P 1,p v p -v p H 1 (Ω) ≤ C 1,p h θ1,p v p H κ+1 (Ω) , ∀v p ∈ H κ+1 (Ω), ∀h ∈ (0, h * );
(H3) ∃C q > 0, ∃θ q ≥ 0 :

- → P q - → v q -- → v q L 2 (Ω) ≤ C q h θq - → v q H κ+1 (div;Ω) , ∀ - → v q ∈ H κ+1 (div; Ω), ∀h ∈ (0, h * ); (H4) ∃C 1,0 > 0, ∃θ 1,0 ≥ 0 : ---→ grad v d p L 2 (Ω) ≤ C 1,0 h -θ1,0 v d p L 2 (Ω) , ∀v d p ∈ H p , ∀h ∈ (0, h * ); (H5) ∃C 0,1 > 0, ∃θ 0,1 ≥ 0 : P p v p -v p H 1 (Ω) ≤ C 0,1 h -θ0,1 P 1,p v p -v p H 1 (Ω) , ∀v p ∈ H 1 (Ω), ∀h ∈ (0, h * ).
Note that in general, (H5) can be deduced from (H4) (see Lemma A.1 in Section A). However this estimate can be strengthened (i.e. θ 0,1 < θ 1,0 ) in many cases in practice, typically with simplicial, regular and quasi-uniform meshes, one has θ 0,1 = 0 thanks to the so-called Aubin-Nitsche trick [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF][START_REF] Sayas | Aubin-Nitsche estimates are equivalent to compact embeddings[END_REF][START_REF] Gatica | A Simple Introduction to the Mixed Finite Element Method: Theory and Applications[END_REF]. It is thus given separately to ensure the optimality of the result.

Basic properties.

It is important to notice the obvious properties between the projectors defined above, and denoted by straight or curly font. These will be useful in the sequel to get from the metric induced by Q to the usual one on L 2 (Ω)

N +1 , and conversely.

• T -1 -→ P q T is a projector from L 2 (Ω) endowed with the inner product

- → v 1 , T - → v 2 L 2 for all - → v 1 , - → v 2 ∈ L 2 (Ω), onto V q ;
• ρP p ρ -1 is a projector from L 2 (Ω) endowed with the inner product

v 1 , ρ -1 v 2 L 2 for all v 1 , v 2 ∈ L 2 (Ω), onto V p ; • T -→ P P P q T -1
is a projector from L 2 (Ω) onto H q ; • ρ -1 P p ρ is a projector from L 2 (Ω) onto H p .

Orthogonality of -→ P q in L 2 (Ω) and P p in L 2 (Ω) imply that ( 18)

- → v - - → P q - → v L 2 ≤ - → v -T -→ P P P q T -1 - → v L 2 , ∀ - → v ∈ L 2 (Ω), v -P p v L 2 ≤ v -ρ -1 P p ρ v L 2 , ∀v ∈ L 2 (Ω), T 1 2 - → v - -→ P P P q - → v L 2 ≤ T 1 2 - → v -T -1 - → P q T - → v L 2 , ∀ - → v ∈ L 2 (Ω), ρ -1 2 (v -P p v) L 2 ≤ ρ -1 2 v -ρP p ρ -1 v L 2 , ∀v ∈ L 2 (Ω).
Thus, for all

- → α q α p ∈ X (19) - → α q α p - -→ P P P q 0 0 P p - → α q α p X ≤ - → α q α p -T -1 - → P q T 0 0 ρP p ρ -1 - → α q α p X ,
and (H1) to (H5) can be written for the curly projectors P p and -→ P P P q in their respective metric, using the upper and lower bounds of the physical parameters ρ and T .

4.2.

From energy variables to co-energy variables. The main results of this section are Theorem 4.2 for E X , and Theorem 4.4 for E H giving the convergence rates in term of those in (H1)-(H5). Theorem 4.4 for the convergence rate of E H emphasizes the interest of the well-known compatibility conditions, appearing e.g. for the finite-dimensional spaces to satisfy an exact sequence, mimicking the de Rham cohomology.

As to avoid introduction of unnecessary notations, the numerical analysis will be carried out making use of the energy variables in X , since this is the natural state space identified in Section 2 [START_REF] Kurula | Linear wave systems on n-D spatial domains[END_REF], thanks to the following (obvious) lemma.

Lemma 4.1. One has

E X (t) := - → α q (t) α p (t) - - → α d q (t) α d p (t) X = T -1 2 0 0 ρ 1 2 - → e (t) e p (t) - - → e d q (t) e d p (t) (L 2 (Ω)) N +1
.

Proof. The proof is straightforward thanks to the constitutive relations -→ α q = T -→ e q and α p = ρ -1 e p .

4.3.

The state absolute error.

Theorem 4.2. Let κ > 0 be an integer and Ω be of class C κ+2 . There exists a constant C X > 0 such that for all T > 0, all initial data

- → α q 0 α p 0 ∈ Z κ , all u ∈ C 2 ([0, ∞); U κ ) such that u(0) = γ ⊥ T - → α q 0 , and all h ∈ (0, h * ) (20) E X (t) ≤ -→ P P P q 0 0 P p - → α q 0 α p 0 - - → α d q (0) α d p (0) X + C X max{1, T } h θ * - → α q α p L ∞ ([0,T ];Zκ) + C X T h -θ1,0 u -u d L ∞ ([0,T ];U ) , ∀t ∈ [0, T ], where (21)
θ * := min {θ 1,p -θ 0,1 ; θ p -θ 1,0 ; θ q -θ 1,0 } .

Proof. For the sake of readability, two technical lemmas are proved in Appendix A.

Remark first that from (H0),

- → α q α p ∈ C([0, ∞); Z κ ),
- → α q α p - - → α d q α d p X ≤ - → α q α p - -→ P P P q 0 0 P p - → α q α p X E1 + -→ P P P q 0 0 P p - → α q α p - - → α d q α d p X E2 .
The strategy of the proof proceeds in five steps.

• In step 1, the convergence of the first term E 1 on the right-hand side of ( 22) is proved thanks to [START_REF] Egger | On Structure-Preserving Model Reduction for Damped Wave Propagation in Transport Networks[END_REF]. • In step 2, Lemma A.2 is applied in order to get the exact value of 1 2

d dt E 2 2 = E 2 d dt E 2 ,
in terms of L 2 (Ω)-inner products and boundary duality bracket.

• In step 3, Cauchy-Schwarz inequality and coarse bounds are used so that

E 2 d dt E 2 ≤ E 3 E 2 .
Dividing by

E 2 > 0 leads to d dt E 2 ≤ E 3 .
• In step 4, Lemma A.3 is applied in order to estimate E 3 .

• Finally, in step 5, the inequality obtained in step 4 is integrated in time, and all the estimates are gathered to conclude.

Step 1 From ( 19) and using the lower bound T -for T and the upper bound ρ + for ρ, one gets (with the definition of the weighted norm on X )

- → α q α p - -→ P P P q 0 0 P p - → α q α p X ≤ - → α q α p -T -1 - → P q T 0 0 ρP p ρ -1 - → α q α p X = T -1 T - → α q ρρ -1 α p -T -1 - → P q T 0 0 ρP p ρ -1 - → α q α p X = T -1 2 - → I q - - → P q 0 0 √ ρ (I p -P p ) T - → α q ρ -1 α p L 2 (Ω)×L 2 (Ω) ≤ 1 √ T- T - → α q - - → P q T - → α q L 2 (Ω) + ρ + ρ -1 α p -P p ρ -1 α p L 2 (Ω) .
Applying (H1) with ρ -1 α p ∈ H κ+1 (Ω) and (H3) with T -→ α q ∈ H κ+1 (div; Ω) leads to

(23) E 1 ≤ max C q T - h θq , ρ + C p h θp - → α q α p Zκ .
Step 2 Applying Lemma A.2 leads to 1 2 d dt

-→ P P P q 0 0 P p - → α q α p - - → α d q α d p 2 X = ---→ grad ρ -1 (α p -P p α p ) , T -→ P P P q - → α q -- → α d q L 2 (Ω) -T - → α q - -→ P P P q - → α q , ---→ grad ρ -1 P p α p -α d p L 2 (Ω) + u -u d , γ 0 ρ -1 P p α p -α d p U ,Y .
Step 3 From Cauchy-Schwarz inequality and the continuity of the Dirichlet trace operator on H 1 (Ω)

1 2 d dt -→ P P P q 0 0 P p - → α q α p - - → α d q α d p 2 X ≤ T 1 2 ---→ grad ρ -1 (α p -P p α p ) L 2 (Ω) T 1 2 -→ P P P q - → α q -- → α d q L 2 (Ω) + T - → α q - -→ P P P q - → α q L 2 (Ω) ---→ grad ρ -1 P p α p -α d p L 2 (Ω) + C D u -u d U ρ -1 P p α p -α d p L 2 (Ω) + C D u -u d U ---→ grad ρ -1 P p α p -α d p L 2 (Ω) . But ρ -1 P p α p -α d p ∈ V p , thus ---→ grad ρ -1 P p α p -α d p L 2 (Ω)
can be estimated by (H4),

---→ grad ρ -1 P p α p -α d p L 2 (Ω) ≤ C 1,0 √ ρ - h -θ1,0 ρ -1 2 P p α p -α d p L 2 (Ω)
, where we have used the lower bound ρ -for ρ. This leads to 1 2 d dt

-→ P P P q 0 0 P p - → α q α p - - → α d q α d p 2 X ≤ T 1 2 ---→ grad ρ -1 (α p -P p α p ) L 2 (Ω) T 1 2 -→ P P P q - → α q -- → α d q L 2 (Ω) + C 1,0 √ ρ - h -θ1,0 T - → α q - -→ P P P q - → α q L 2 (Ω) ρ -1 2 P p α p -α d p L 2 (Ω) + C D √ ρ - u -u d U ρ -1 2 P p α p -α d p L 2 (Ω) + C D C 1,0 √ ρ - h -θ1,0 u -u d U ρ -1 2 P p α p -α d p L 2 (Ω)
.

Gathering T 1 2

-→ P P P q - → α q -- → α d q L 2 (Ω)
and ρ -1 2

P p α p -α d p L 2 (Ω)
gives the desired X -norm:

1 2 d dt -→ P P P q 0 0 P p - → α q α p - - → α d q α d p 2 X ≤ T 1 2 ---→ grad ρ -1 (α p -P p α p ) L 2 (Ω) + C 1,0 √ ρ - h -θ1,0 T - → α q - -→ P P P q - → α q L 2 (Ω) + C D √ ρ - 1 + C 1,0 h -θ1,0 u -u d U -→ P P P q 0 0 P p - → α q α p - - → α d q α d p X
.

Dividing both sides by

-→ P P P q 0 0 P p - → α q α p - - → α d q α d p X
, we finally get

d dt E 2 ≤ E 3 , with E 3 := T 1 2 ---→ grad ρ -1 (α p -P p α p ) L 2 (Ω) + C 1,0 √ ρ - h -θ1,0 T - → α q - -→ P P P q - → α q L 2 (Ω) + C D √ ρ - 1 + C 1,0 h -θ1,0 u -u d U .
Step 4 Using the upper bound T + for T , we get

E 3 := √ T + ---→ grad ρ -1 (α p -P p α p ) L 2 (Ω) + C 1,0 √ ρ - h -θ1,0 T - → α q - -→ P P P q - → α q L 2 (Ω) + C D √ ρ - 1 + C 1,0 h -θ1,0 u -u d U .
From Lemma A.3, one has

(24) E 3 ≤ √ T + C 0,1 C 1,p h θ1,p-θ0,1 + ρ + C 1,0 C p √ ρ - h θp-θ1,0 ρ -1 α p H κ+1 (Ω) + C 1,0 √ ρ - h -θ1,0 T - → α q - -→ P P P q - → α q L 2 (Ω) + C D √ ρ - 1 + C 1,0 h -θ1,0 u -u d U .
It remains to estimate T -→ α q --→ P P P q -→ α q L 2 (Ω)

. Thanks to the upper bound T + for T and the third line of (18), we have

T - → α q - -→ P P P q - → α q L 2 (Ω) ≤ √ T + T 1 2 - → α q -T -1 - → P q T - → α q L 2 (Ω)
, or in other words

T - → α q - -→ P P P q - → α q L 2 (Ω) ≤ √ T + T -1 2 T - → α q - - → P q T - → α q L 2 (Ω)
.

With the lower bound T -for T and (H3) with -→ v q = T -→ α q , this leads to

T - → α q - -→ P P P q - → α q L 2 (Ω) ≤ √ T + T - C q h θq T - → α q H κ+1 (div;Ω)
.

By injecting the latter estimate into (24), we get

E 3 ≤ √ T + C 0,1 C 1,p h θ1,p-θ0,1 + ρ + C 1,0 C p √ ρ - h θp-θ1,0 ρ -1 α p H κ+1 (Ω) + √ T + C 1,0 C q T -ρ - h θq-θ1,0 T - → α q H κ+1 (div;Ω) + C D √ ρ - 1 + C 1,0 h -θ1,0 u -u d U ,
which gives, with a rough majoration, the existence of a constant C 3 > 0 such that

E 3 ≤ C 3 h min{θ1,p-θ0,1 ; θp-θ1,0 ; θq-θ1,0} - → α q α p Zκ + C 3 h -θ1,0 u -u d U ,
for all h small enough.

Step 5 By integrating d dt E 2 ≤ E 3 between 0 and T , the latter inequality gives 23) and ( 25) into [START_REF] Farle | A port-Hamiltonian finite-element formulation for the Maxwell equations[END_REF], noticing that θ p ≥ θ p -θ 1,0 and θ q ≥ θ q -θ 1,0 , gives the desired result for all h small enough.

(25) E 2 (t) ≤ E 2 (0) + C 3 T h min{θ1,p-θ0,1 ; θp-θ1,0 ; θq-θ1,0} - → α q α p L ∞ ([0,T ];Zκ) + C 3 T h -θ1,0 u -u d L ∞ ([0,T ];U ) . Substituting (

The Hamiltonian error.

In this subsection, the numerical analysis focuses on the Hamiltonian, main object of interest in the port-Hamiltonian framework.

The next corollary follows easily from Theorem 4.2, despite it does not give the expected optimal convergence rate: twice that of the state space absolute error.

(26) E H (t) ≤ - → α q α p L ∞ ([0,T ];X ) + E X (t) 2 E X (t), ∀t ∈ [0, T ].
Proof. It is straightforward that

E H = 1 2 - → α q α p + - → α d q α d p , - → α q α p - - → α d q α d p X ≤ 1 2 - → α q α p + - → α d q α d p X E X . But - → α q (t) α p (t) + - → α d q (t) α d p (t) X ≤ 2 - → α q α p L ∞ ([0,T ];X ) + E X (t), ∀t ∈ [0, T ],
which ends the proof of (26).

In the following theorem, it is proved that compatibility conditions between H q , H p and H ∂ , including but not restricted to those leading to the preservation of the de Rham cohomology (see e.g. the commutative diagrams [10, page 116]) at the discrete level, lead to a far better result for E H . Theorem 4.4. Under the assumptions of Theorem 4.2, assume furthermore that

• ---→ grad vp-Ppvp ρ , - → v d q L 2 (Ω) = 0, for all - → v d q ∈ H q , v p ∈ H 1 (Ω); • T - → v q - -→ P P P q - → v q , ---→ grad v d p L 2 (Ω) = 0, for all v d p ∈ H p , - → v q ∈ L 2 (Ω); • u -u d , γ 0 v d p L 2 (∂Ω) = 0, for all u ∈ L 2 (∂Ω), u d approximation of u in H ∂ , v d p ∈ H p ; • u d , γ 0 v d p -Ppvp ρ L 2 (∂Ω) = 0, for all u d ∈ H ∂ , v p ∈ H 1 (Ω), v d p approxi- mation of v p in H p . Then E H (t) -E H (0) = 1 2 E X (t) 2 -E X (0) 2 , ∀t ∈ [0, T ]. Proof. Obviously E X 2 = - → α q α p 2 X -2 - → α q α p , - → α d q α d p X + - → α d q α d p 2 X
. Hence ( 27)

E H = 1 2 E X 2 + - → α q α p , - → α d q α d p X - - → α d q α d p 2 X .
From ( 10) and ( 11) together with Lemma 4.1, it is straightforward that d dt

- → α q α p , - → α d q α d p X = ---→ grad ρ -1 α p - P p α p ρ , T - → α d q L 2 (Ω) + T -→ P P P q - → α q -T - → α q , ---→ grad α d p ρ L 2 (Ω) + u, γ 0 α d p ρ L 2 (∂Ω) + u d , γ 0 P p α p ρ L 2 (∂Ω)
, which becomes, thanks to the assumptions on the q-and p-type families (remember

that - → α d q ∈ V q := T -1 H q and α d p ∈ V p := ρH p ) d dt - → α q α p , - → α d q α d p X = u, γ 0 α d p ρ L 2 (∂Ω) + u d , γ 0 P p α p ρ L 2 (∂Ω)
.

Since by Proposition 3.1

d dt - → α d q α d p 2 X = 2 u d , y d L 2 (∂Ω) ,
it can be deduced that d dt

- → α q α p , - → α d q α d p X - - → α d q α d p 2 X = u -u d , γ 0 α d p ρ L 2 (∂Ω) + u d , γ 0 P p α p -α d p ρ L 2 (∂Ω)
. Now, thanks to the assumptions involving the boundary finite element families, an integration in time from 0 to t gives

- → α q α p , - → α d q α d p X - - → α d q α d p 2 X = - → α q 0 α p 0 , - → α d q (0) α d p (0) X - - → α d q (0) α d p (0) 2 X ,
and the result follows from [START_REF] Gerritsma | Discrete Lie derivative[END_REF] by subtracting E H (0) from both side.

Remark 4.3. Note that the anisotropy and heterogeneity have a non-negligible influence for the validity of Theorem 4.4, as these induce curly projectors -→ P P P q and P p . These modifications of the compatibility conditions for this more general case would be hidden if the analysis was carried out with constant parameters. 4.5. Other causality: the Dirichlet boundary control. Theorem 4.2 has its counterpart for the other causality, already discussed in Section 3.2, where u and y have been switched (notation S). For the sake of space saving, since the proof is quite similar, only the result for the general framework are briefly stated below.

Assuming an existence and regularity result such as (H0) for the case (1)-(2S) (though not covered in [START_REF] Kurula | Linear wave systems on n-D spatial domains[END_REF]), one can easily adapt the proof of Theorem 4.2.

If

• H q is H(div; Ω)-conforming (instead of L 2 (Ω)-conforming); • H p is L 2 (Ω)-conforming (instead of H 1 (Ω)-conforming); • H ∂ is H 1 2 (∂Ω)-conforming, satisfying (H1S) ∃C p > 0, ∃θ p ≥ 0 : P p v p -v p L 2 (Ω) ≤ C p h θp v p H κ+1 (Ω) , ∀v p ∈ H κ+1 (Ω), ∀h ∈ (0, h * );
(H2S) ∃C div,q > 0, ∃θ div,q ≥ 0 :

- → P div,q - → v q -- → v q H(div;Ω) ≤ C div,q h θ div,q - → v q H κ+1 (div;Ω) , ∀ - → v q ∈ H κ+1 (div; Ω), ∀h ∈ (0, h * );
where -→ P div,q is the H(div; Ω)-orthogonal projector from H(div; Ω) onto V q ;

(H3S) ∃C q > 0, ∃θ q ≥ 0 :

- → P q - → v q -- → v q L 2 (Ω) ≤ C q h θq - → v q H κ+1 (div;Ω) , ∀ - → v q ∈ H κ+1 (div; Ω), ∀h ∈ (0, h * ); (H4S) ∃C div,0 > 0, ∃θ div,0 ≥ 0 : div - → v d q L 2 (Ω) ≤ C div→0 h -θ div,0 - → v d q L 2 (Ω) , ∀ - → v d
q ∈ H q , ∀h ∈ (0, h * );

(H5S) ∃C 0,div > 0, ∃θ 0,div ≥ 0 :

- → P q - → v q -- → v q H(div;Ω) ≤ C 0,div h -θ 0,div - → P div,q - → v q -- → v q H(div;Ω) , ∀ - → v q ∈ H(div; Ω), ∀h ∈ (0, h * );
we have the following theorem.

Theorem 4.5. Let κ > 0 be an integer and Ω be of class C κ+2 . There exists a constant C X > 0 such that for all T > 0, all initial data

- → α q 0 α p 0 ∈ Z κ , all u ∈ C 2 ([0, ∞); U κ ) such that u(0) = γ ⊥ T - → α q 0 , and all h ∈ (0, h * ) E X (t) ≤ -→ P P P q 0 0 P p - → α q 0 α p 0 - - → α d q (0) α d p (0) X + C X max{1, T } h θ * - → α q α p L ∞ ([0,T ];Zκ) + C X T h -θ div,0 u -u d L ∞ ([0,T ];U ) , ∀t ∈ [0, T ],
where θ * := min {θ div,q -θ 0,div ; θ q -θ div,0 ; θ q -θ div,0 } . Remark 4.4. A counterpart of Theorem 4.4 should also be possible to prove under suitable additional compatibility conditions.

Optimal orders of conforming finite elements

The purpose of this section is to provide optimal combinations of usual finite elements minimizing the number of degrees of freedom for a given convergence rate, illustrating the abstract estimates obtained in Theorems 4.2 and 4.4. The errors to be analysed when the mesh size parameter h > 0 tends towards 0 are those of the previous section, and will be numerically investigated in the next section.

5.1. Mesh assumptions. These classical assumptions in numerical analysis for usual finite elements (see e.g. [26, pp. 61, 88, 96] and [10, p. 71]) are recalled for the sake of completeness.

The mesh family (T h ) h∈(0,h * ) of Ω will be supposed to be a collection of simplicial, regular and quasi-uniform triangularization of Ω, meaning that (1) It is given by a collection of triangles or tetrahedra, denoted K in the sequel.

(

) If h K > 0 denotes the diameter of K, i.e. h K := max-→ x , - → y ∈K - → x -- → y , and 2
d K is the diameter of the inscribed circle or sphere in K, there exists a constant C > 0, independent of h such that

h K d K ≤ C, ∀K ∈ T h , ∀h ∈ (0, h * ).
The mesh parameter is then defined as h := max K∈T h h K . (3) There exists a constant c > 0 independent of h such that min K∈T h h K ≥ ch for all h ∈ (0, h * ).

Lagrange element.

The three finite element families must be L 2 (Ω)-, H 1 (Ω)-and H 1 2 (∂Ω)-conforming respectively, in order to apply Theorem 4.2. A first easy choice is to take the usual continuous Galerkin finite elements for families q and p, i.e.

H q := - → v d q ∈ C(Ω) N | - → v d q K ∈ (P k (K)) N , ∀K ∈ T h ,
when k ≥ 1, or the piecewise constant functions when k = 0.

H p := v d p ∈ C(Ω) | v d p K ∈ P (K), ∀K ∈ T h ,
where ≥ 1.

In the above definition, P j (K) is the Lagrange finite element of order j made of all polynomials of degree less or equal to j on K.

The space H p is known as continuous Galerkin of order as we impose continuity of basis functions. The space H q is the vectorial counterpart of H p , at order k. In the sequel, we will refer to these spaces via CG and CG k respectively.

For the discretization space at the boundary H ∂ , discontinuous Galerkin finite elements of order m ≥ 0, denoted DG m , are chosen

H ∂ := v d ∂ ∈ L ∞ (∂Ω) | v d ∂ E ∈ P m (E), ∀E,
edges of K ∈ T h located at the boundary .

They are indeed in H 1 2 (∂Ω) for all m ≥ 0. Practically, it consists in taking the Dirichlet trace of discontinuous Galerkin of order m on the whole domain, i.e. keeping only degrees of freedom located at the boundary.

All error estimates are well-known and can be found e.g. in [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF][START_REF] Gatica | A Simple Introduction to the Mixed Finite Element Method: Theory and Applications[END_REF] and references therein for the global interpolation operators. Obviously, these estimates hold true for the orthogonal projectors. Then (H1)-(H2)-( H3)-(H4)-(H5) read (H1L) ∃C p > 0,

P p v p -v p L 2 (Ω) ≤ C p h +1 v p H +1 (Ω) , ∀v p ∈ H +1 (Ω), ∀h ∈ (0, h * ); (H2L) ∃C 1,p > 0, P 1,p v p -v p H 1 (Ω) ≤ C 1,p h v p H +1 (Ω) , ∀v p ∈ H +1 (Ω), ∀h ∈ (0, h * ).
The following estimate requires more attention. Nevertheless, a careful analysis using (H1L), [START_REF] Foiaş | Remarques sur les équations de Navier-Stokes stationnaires et les phénomènes successifs de bifurcation[END_REF]Proposition 1.4], and a density argument, leads to (H3L) ∃C q > 0,

- → P q - → v q -- → v q L 2 (Ω) ≤ C q h k+1 - → v q H k+1 (div;Ω) , ∀ - → v q ∈ H k+1 (div; Ω), ∀h ∈ (0, h * ). (H4L) ∃C 1,0 > 0, ---→ grad v d p L 2 (Ω) ≤ C 1,0 h -1 v d p L 2 (Ω) , ∀v d p ∈ H p , ∀h ∈ (0, h * ),
and finally (H5L) ∃C 0,1 > 0,

P p v p -v p H 1 (Ω) ≤ C 0,1 P 1,p v p -v p H 1 (Ω) , ∀v p ∈ H 1 (Ω), ∀h ∈ (0, h * ).
Remark that, with the choice of H ∂ , a straightforward but tedious exercise, using lifting operators, Bramble-Hilbert Theorem [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF], continuity of trace operators and quasi-uniform hypothesis give that [START_REF] Gerritsma | Mimetic Spectral Element Method for Anisotropic Diffusion[END_REF] u

-u d H -1 2 (∂Ω) ≤ C u h m+1 u H m-1 2 (∂Ω)
, ∀u ∈ H m-1 2 (∂Ω), at the boundary.

The following holds true.

Theorem 5.1. Let κ > 0 be an integer, T > 0 a final time,

- → α q 0 α p 0 ∈ Z κ the initial data, u ∈ C 2 ([0, ∞); H κ-1 2 (∂Ω)
) the control, and

- → α d q (0) α d p (0)
and u d their respective continuous and discontinuous Galerkin interpolations in V q × V p and H ∂ given by the finite elements

(CG k ) N × CG × DG m .
There exists a constant C > 0, independent of T > 0, -→ α q 0 α p 0 , and u, such that for all h small enough and all t ∈ [0, T ]

(29) E X (t) ≤ C max{1, T } h min{ ; k ; m} × - → α q α p L ∞ ([0,T ];Zκ) + u L ∞ ([0,T ];H κ-1 2 (∂Ω))
.

Furthermore, the optimal order is κ, obtained with k = κ, = κ and m = κ -1.

Proof. Since κ ≥ 1, H κ-1 2 (∂Ω) = U κ ⊂ H 1 2 (
∂Ω), and u can indeed be approximated in H ∂ .

From (H0), the solution to ( 1)-( 2) belongs to Z κ continuously in time. Recall that this means

- → α q ∈ T -1 H κ+1 (div; Ω) := - → v ∈ L 2 (Ω) | T - → v ∈ H κ (Ω), div T - → v ∈ H κ (Ω) ,
and α p ∈ ρH κ+1 (Ω). Hence, following (H1L)-( H2L)-( H3L)-( 28), the order of the finite element families satisfies

θ p = + 1 ≤ κ + 1, θ 1,p = ≤ κ, θ q = k + 1 ≤ κ + 1, θ u := m + 1 ≤ κ.
From (H4L)-(H5L)

θ 1,0 = 1, θ 0,1 = 0.
One deduces the convergence rate of E X thanks to ( 20) and ( 21), i.e. it is given by

θ * = min { ; k ; m + 1} ,
where we have used (H1L) and (H3L) for the approximation of the initial data. Now, taking into account the maximal regularities given by (H0) (and the assumed regularity on u) leads to the maximal rate min {κ ; κ ; κ} = κ. Finally, one gets this maximal order with the minimal number of degrees of freedom when we take = k = m + 1 = κ.

5.3.

Other finite element families. Following [10, Proposition 2.5.4.], estimate (H3), and thus Theorem 4.2, hold true for many usual H(div; Ω)-conforming families (hence L 2 (Ω)-conforming as required, or even with curl-conforming finite element), namely: Raviart-Thomas RT k (for an introduction to this important class of finite element, see e.g. [START_REF] Gatica | A Simple Introduction to the Mixed Finite Element Method: Theory and Applications[END_REF]), Brezzi-Douglas-Marini BDM k and discontinuous Galerkin finite elements DG k . Proposition 5.2. Let κ > 0 be an integer, T > 0 a final time,

- → α q 0 α p 0 ∈ Z κ the initial data, u ∈ C 2 ([0, ∞); H κ-1 2 (∂Ω)
) the control, and

- → α d q (0) α d p (0)
and u d their respective interpolations. The optimal rate of convergence is reached with V q ×V p × H ∂ given by

DG κ-1 × CG κ × DG κ-1 , RT κ × CG κ × DG κ-1 , BDM κ × CG κ × DG κ-1 , CG κ × CG κ × CG κ , DG κ-1 × CG κ × CG κ , RT κ × CG κ × CG κ , BDM κ × CG κ × CG κ ,
all of them leading to the same convergence rate κ.

Proof. This is a direct application of Theorem 4.2.

Remark 5.1. Care must be taken with the subscript of RT element, which may differ from one source to another, depending on how the lowest order is denoted: either RT 0 or RT 1 . In this paper, we stick to the definition given in FEniCS [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF], the software being used for the simulations ran in Section 6, and denote the lowest order by RT 1 .

Numerical study of the convergence rate in 2D

In this section, simulations are performed to illustrate our results. More precisely, we intend to verify if the convergence rates are indeed those proved in Theorem 5.1 and claimed in Section 5.3.

6.

1. An analytical solution. In order to study the convergence rate, we propose to focus on a 2D toy model, isotropic and heterogeneous, for which an analytical solution is known. This choice is made to avoid the computation of a reference solution.

Let us consider Ω = (0, 1)×(0, 1). The physical parameters are ρ ≡ 1 and

T ≡ I. Denoting f (t) := 2 sin √ 2t + 3 cos √ 2t , we define - → α q := f (t) -sin(x) sin(y) cos(x) cos(y) , α p := d dt f (t) cos(x) sin(y), ∀(x, y) ∈ Ω, t ≥ 0, and
u(t) :=        -f (t) cos(x), ∀(x, y) ∈ (0, 1) × {0}, -f (t) sin(1) sin(y), ∀(x, y) ∈ {1} × (0, 1), f (t) cos(x) cos(1), ∀(x, y) ∈ (0, 1) × {1}, 0, ∀(x, y) ∈ {0} × (0, 1).
Then,

- → α q α p is a C ∞ ([0, ∞); C ∞ (Ω)
)-solution to the wave equation written as a pHs ( 6)-( 7). The choice of sine and cosine functions has been made to avoid exact interpolation in the polynomial finite element spaces of high order.

The Hamiltonian is easily obtained for all t ≥ 0

H - → α q (t), α p (t) = 1 8 d dt f (t) 2 1 -(sin(1) cos(1)) 2 + 1 8 (f (t)) 2 × (1 + sin(1) cos(1)) 2 + (1 -sin(1) cos(1)) 2 .
Figure 1. Convergence of E X as a function of h at order 1, with the optimal choice of finite element families. See values in bold on Table 1.

6.2. Simulations. In this section, the following procedure is proposed: H q , H p and H ∂ are varying according to many ranges of finite element families, and all combinations are tested to analyse the behaviour of the convergence rate. The simulations are performed using FEniCS [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF], with a Crank-Nicolson scheme in time t ∈ (0, 0.5) [START_REF] Abhyankar | PETSc/TS: A modern scalable ODE/DAE solver library[END_REF]. The time step is chosen small enough according to κ, in order to ensure that the error is driven by the mesh size h and not the time 3.

step 1 dt = 10 -5 . Note that Crank-Nicolson scheme being an implicit time-stepping scheme, it requires the resolution of a linear system at each time step: difficulties may arise due to the behaviour of condition numbers of the matrices involved, when h → 0, which is often the case for finite element methods. Since the purpose of 1 This is not a kind of CFL condition due to the PFEM, but a matter of space discretization analysis. We only have access to the total error O(h κ ) + O(∆t 2 ). Taking a sufficiently small time step allows for the analysis of this error as a function of h, to be able to observe the order κ. Note in particular that Crank-Nicolson scheme is unconditionally stable and symplectic for linear systems.

Table 1.

Convergence at order 1 for E X obtained for different combinations of finite element families. The columns correspond to the q-type vector-valued variables, and the rows to the p-type scalar ones. The orders in bold are the optimal orders given by Theorem 5.1 and Proposition 5.2, which give the expected convergence rate with the lowest complexity. obtained for different combinations of finite element families. The columns correspond to the q-type vector-valued variables, and the rows to the p-type scalar ones. The orders in bold are the optimal orders given by Theorem 5.1 and Proposition 5.2, which give the expected convergence rate with the lowest complexity. CG1 boundary finite elements on ∂Ω Table 3. Convergence at order 3 for E X obtained for different combinations of finite element families. The columns correspond to the q-type vector-valued variables, and the rows to the p-type scalar ones. The orders in bold are the optimal orders given by Theorem 5.1 and Proposition 5.2, which give the expected convergence rate with the lowest complexity. this work is semi -discretization in space, we do not go deeper into that direction and refer the interested reader to e.g. [START_REF] Fried | Bounds on the extremal eigenvalues of the finite element stiffness and mass matrices and their spectral condition number[END_REF][START_REF] Ern | Evaluation of the condition number in linear systems arising in finite element approximations[END_REF] for that specific problem, and to [START_REF] Kotyczka | Discrete-time port-Hamiltonian systems: A definition based on symplectic integration[END_REF] for discrete-time finite-dimensional port-Hamiltonian systems. The tests (namely 447 combinations, with 8 different values for h) have been run on a personal computer (Intel Core I7 processor, 24GB of RAM). All the convergence rates are presented on Tables 1-2-3 for the absolute error E X . The associated figures show the convergence rates proven in Theorem 4.2 and Proposition 5.2 for the optimal choices of finite element families among all our tests.

DG

In Theorem 4.2, H p is assumed to be H 1 (Ω)-conforming. Looking at Tables 1-2-3, this seems indeed necessary to ensure convergence. In every test cases, the rate is negative when H p is given by discontinuous Galerkin finite elements, no matter the order (we gather lines for DG , = 0, 1, 2, 3 since they give the same results). It has to be noted that assumption (H2) is not satisfied in these cases: Theorem 4.2 does not apply.

6.3. About the convergence rate of the Hamiltonian error E H . So far, compatibility conditions have not been taken into account between H q , H p and H ∂ . In other words, only conforming assumptions have been made, and optimal rates have then been deduced. However, pHs are strongly structured, and in particular, the de Rham cohomology should be respected to improve the efficiency of the PFEM. As a motivation, it is remarkable on Tables 4-5-6 that E H does not converge at the same rate as E X , as stated in Corollary 4.3, but at twice its order in a various number of cases.

On Table 4, one can see that the Hamiltonian convergence rates when boundary functions are approximated by Discontinuous Galerkin finite elements of order 0, DG 0 , never reach order 2, i.e. twice the optimal convergence rate of the state error according to Theorem 5.1. This leads us to conclude that compatibility conditions of Theorem 4.4 are never met for the combinations of finite elements presented on Table 4.

On the contrary, increasing by one order the approximation for boundary terms, both for Discontinuous and Continuous Galerkin DG 1 and CG 1 , leads to an order 2 in most cases, as seen on Table 5. Those corresponding to a convergence rate of order 1 for the state error in Table 2 might indeed satisfy the compatibility condition of Theorem 4.4.

Analogously, increasing again the approximation at the boundary, i.e. taking DG 2 or CG 2 finite elements, allows for a Hamiltonian convergence rate reaching order 4, inviting us to conjecture that compatibility conditions of Theorem 4.4 are met for boldface convergence rate of Table 6. Remark that we do not have a sufficiently small time step to be able to numerically observe order 6, if any.

To conclude, as pHs deal with a Hamiltonian functional, which can be seen as the primary object, the structure-preserving discretization should mean that H has to be accurately discretized for both the power balance (PFEM) and the value of H in R (compatibility conditions). Together, the PFEM and the compatibility conditions seem to achieve this, making use of the Finite Element Method only. They give the maximal precision with the minimal number of degrees of freedom. Furthermore, the number of degrees of freedom at the boundary being very low in comparison to those of q-and p-type, it clearly appears that the choice of boundary finite elements proves crucial to reach the expected convergence rates with respect to the finite elements chosen within the domain Ω (as an example, compare the rate obtained with the discretization CG 2 × CG 3 for -→ α q and α p on Table 5 with the same on Table 6). An analytical solution is then given by: w(t, x) = cos(3t -x + 2y). More precisely

DG 3 CG 1 CG 2 CG 3 BDM 1 BDM 2 BDM 3 RT 1 RT 2 RT 3 DG -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 CG 1 1
DG 3 CG 1 CG 2 CG 3 BDM 1 BDM 2 BDM 3 RT 1 RT 2 RT 3 DG -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 CG 1 1
- → α q = -1 2 sin(3t -x + 2y), α p = 3 sin(3t -x + 2y).
Convergence rates for optimal combinations of finite elements are given on Figure 5. For instance in the port-Hamiltonian formalism, at the continuous level, the dissipation does not appear as a positive bounded operator R in the dynamics in general. This means that the dynamics of the port-Hamiltonian systems is not necessarily governed by an operator of the form (J -R)Q (even with a lot of work involving lifting operators to fit this formulation, R will not be a bounded operator). In other words, the dissipativity is hidden in the domain of the unbounded operator J . However, it is expected to recover a finitedimensional pHs driven by a matrix of the form J d -R d Q d at the discrete level. PFEM provides a very easy and convenient way to describe this positive matrix R d when dealing with admittance or impedance boundary conditions, i.e. absorbing boundary condition in the PDE terminology. The strategy relies on the use of a suitable output feedback law on the finite-dimensional pHs obtained by PFEM in the previous sections (see [START_REF] Serhani | Partitioned Finite Element Method for port-Hamiltonian systems with Boundary Damping: Anisotropic Heterogeneous 2-D wave equations[END_REF][START_REF] Serhani | A Partitioned Finite Element Method for the Structure-Preserving Discretization of Damped Infinite-Dimensional Port-Hamiltonian Systems with Boundary Control[END_REF] where this original strategy has been first proposed and explained in full details). Furthermore, this construction gives a well-understood structure to the matrix R d , which turns out to be of low rank: at most the dimension of V ∂ as expected, since the damping is applied at the boundary only.

6.5.1. Discretization. More precisely, the following boundary condition is considered, instead of (2), for the admittance boundary condition

(2Y)      v(t, - → x ) = Y ( - → x )∂ t w(t, - → x ) + T ( -→ x ) ---→ grad(w(t, - → x )) - → n ( - → x ), ∀ - → x ∈ ∂Ω, t ≥ 0, y(t, - → x ) = ∂ t w(t, - → x ), ∀ - → x ∈ ∂Ω, t ≥ 0.
and the following one instead of (2), for the impedance boundary condition

(2Z)        v(t, - → x ) = ∂ t w(t, - → x ) + Z( - → x ) T ( - → x ) ---→ grad(w(t, - → x )) - → n ( - → x ), ∀ - → x ∈ ∂Ω, t ≥ 0, y(t, - → x ) = T ( - → x ) ---→ grad(w(t, - → x )) - → n ( - → x ), ∀ - → x ∈ ∂Ω, t ≥ 0.
where both the admittance Y and the impedance Z are positive and belong to L ∞ (∂Ω) and v or v are the external inputs.

Remark 6.1. It is clear that (2Y) and (2Z) generalize ( 2) and (2S) respectively. Nevertheless, as mentioned above, PFEM does not apply straightforwardly in these more general cases (think about the use of Green's formula (5) at the beginning of the strategy). The proposed alternative to construct the finite-dimensional dissipative pHs using an output feedback laws seems an elegant way to achieve our goal.

It is easy to write the following relations: between u, y and v (the new control) from ( 2) and (2Y)

u(t, - → x ) = v(t, - → x ) -Y ( - → x )y(t, - → x ), ∀ - → x ∈ ∂Ω, t ≥ 0,
or between u, y and v (the new control) from (2S) and (2Z)

(30) u(t, - → x ) = v(t, - → x ) -Z( - → x ) y(t, - → x ), ∀ - → x ∈ ∂Ω, t ≥ 0.
Using a weak formulation, these equalities read in matrix form

M ∂ u(t) = M ∂ v(t) -Y y(t), ∀t ≥ 0, or M ∂ u(t) = M ∂ v(t) -Z y(t), ∀t ≥ 0,
respectively, where

Y := ∂Ω Y ( - → s )Ψ( - → s ) (Ψ( - → s )) d - → s , or Z := ∂Ω Z( - → s )Ψ( - → s ) (Ψ( - → s )) d - → s .
As presented in [50, Remark 2.], this procedure indeed gives rise to finite-dimensional Dirac structures, introducing extra resistive ports, and leading to a pHDAE.

Simulation results.

As a worked-out example, let us consider a fully heterogeneous and anisotropic case, with boundary control and boundary damping. The aim is to illustrate how the structure-preserving scheme can be appreciate on the Hamiltonian behaviour and the different kind of energies present in the system (potential, kinetic, supplied and damped).

Let us consider each part of the energy and their sum. The preservation of the physical meaning supposes that the exchanges of energy (e.g. potential to kinetic and vice-versa, boundary-supplied/taken energy to the system and damped into internal energy) must result in the preservation of the first principle of thermodynamics: the sum of all energies must be constant over time. More precisely, let us define the potential energy and the first principle of thermodynamics implies that E(t) ≡ E(0) for all t ≥ 0. The simulation is performed on the time interval (0, 3) with a time step dt = 10 -4 . The spatial discretization is RT 1 × CG 1 × CG 1 . The time solver is Assimulo (IDA SUNDIALS) [START_REF] Andersson | Assimulo: A unified framework for ODE solvers[END_REF]. The integration in time to compute S and D are done using the midpoint rule. We can appreciate on Figure 6 how the total energy E remains constant over time, as physically expected, thanks to the PFEM.

E Pot (t) := 1 2 Ω - → α q (t, - → x ) T (x, y) - → α q (t, - → x ) d - → x

Conclusion and perspectives

In this work, the numerical analysis of the Partitioned Finite Element Method applied to the anisotropic and heterogeneous boundary-controlled-and-observed Ndimensional wave equation has been carried out. This recent structure-preserving method [START_REF] Cardoso-Ribeiro | A structure-preserving Partitioned Finite Element Method for the 2D wave equation[END_REF][START_REF] Cardoso-Ribeiro | A partitioned finite element method for power-preserving discretization of open systems of conservation laws[END_REF] allows the direct construction of a finite-dimensional port-Hamiltonian system, the underlying Dirac structure of which mimicks the infinite-dimensional Stokes-Dirac structure. This property allows a very accurate discretization of the power balance satisfied by the system. Furthermore, it has been shown that under compatibility conditions resembling those allowing for discrete de Rham complexes, the discrete Hamiltonian converges very fastly toward the continuous one, strengthening the interest of the PFEM, since the versatility of port-Hamiltonian systems aims precisely at modelling the exchange of energies between sub-systems. As an illustration of our main theorems, we have performed 2D simulations on a case where an analytical solution is known, with a non-homogeneous boundary condition: the boundary control, for both convex and non-convex domain. A wide range of usual finite element have been tested. Moreover, impedance-like boundary damping have been carried out to illustrate the structure-preserving property of the PFEM [START_REF] Serhani | Partitioned Finite Element Method for port-Hamiltonian systems with Boundary Damping: Anisotropic Heterogeneous 2-D wave equations[END_REF][START_REF] Serhani | A Partitioned Finite Element Method for the Structure-Preserving Discretization of Damped Infinite-Dimensional Port-Hamiltonian Systems with Boundary Control[END_REF].

Several questions remain open, the first one being the case of mixed boundary conditions. Two approaches have been proposed in [START_REF] Brugnoli | Partitioned Finite Element Method for Power-Preserving Structured Discretization with Mixed Boundary Conditions[END_REF]: a domain decomposition followed by a gyrator interconnection between the two sub-systems, and the use of Lagrange multipliers. The former has the great advantage of remaining an ODE, while the latter transforms into a Differential Algebraic Equation (DAE). The numerical analysis of both alternatives will require deeper investigation.

A second interesting question is the problem of symplectic integration of pH-DAE, which naturally arises in various situations such as the aforementioned mixed boundary conditions, or for the heat equation [START_REF] Serhani | Anisotropic heterogeneous n-D heat equation with boundary control and observation: I. Modeling as port-Hamiltonian system[END_REF][START_REF] Serhani | Anisotropic heterogeneous n-D heat equation with boundary control and observation: II. Structure-preserving discretization[END_REF]. A scheme has been recently proposed in [START_REF] Mehrmann | Structure-Preserving Discretization for Port-Hamiltonian Descriptor Systems[END_REF] for this purpose.

Finally, in assumption (H0) allowing for regular solutions, X κ is assumed to be a Hilbert space. It is clear that this requires some regularity assumptions on the physical parameters. An interesting future work would be to design a structurepreserving scheme when the physical parameters ρ and T are no more regular, e.g. using mollifying methods: refining the constants in the estimates as function of those parameters, as done in this work, seems to be a first step in order to tackle this problem.

Thanks to the orthogonality of -→ P P P q 0 0 P p in X , we have

∂ t - → α q -∂ t - → α d q , T
-→ P P P q -→ α q --→ α d q L 2 (Ω)

+ ∂ t α p -∂ t α d p , ρ -1 P p α p -α d p L 2 (Ω)
= ∂ t -→ P P P q -→ α q -∂ t -→ α d q , T -→ P P P q -→ α q --→ α d q L 2 (Ω) ≤ C 0,1 C 1,p h θ1,p-θ0,1 + ρ + C 1,0 C p √ ρ - h θp-θ1,0 ρ -1 α p H κ+1 (Ω) .

+
Proof. Writing ---→ grad ρ -1 (α p -P p α p ) , the first term on the right-hand side is bounded thanks to (H5), with v p = ρ -1 α p ∈ H 1 (Ω), ---→ grad ρ -1 (α p -P p α p )

L 2 (Ω) ≤ C 0,1 h -θ0,1 P 1,p ρ -1 α p -ρ -1 α p H 1 (Ω) + ---→ grad P p ρ -1 α p -ρ -1 P p α p L 2 (Ω)
, and by (H2), still with v p = ρ -1 α p , ---→ grad ρ -1 (α p -P p α p )

L 2 (Ω)

≤ C 0,1 C 1,p h θ1,p-θ0,1 ρ -1 α p H κ+1 (Ω) + ---→ grad P p ρ -1 α p -ρ -1 P p α p L 2 (Ω)

.

Since (P p ρ -1 α p -ρ -1 P p α p ) ∈ V p , hypothesis (H4) gives ---→ grad ρ -1 (α p -P p α p )

L 2 (Ω)

≤ C 0,1 C 1,p h θ1,p-θ0,1 ρ -1 α p H κ+1 (Ω) +C 1,0 h -θ1,0 P p ρ -1 α p -ρ -1 P p α p L 2 (Ω) .

Let us focus now on P p ρ -1 α p -ρ -1 P p α p L 2 (Ω) to conclude. Since ρ -1 P p ρ is a projector from L 2 (Ω) onto V p , one has P p ρ -1 α p = ρ -1 P p ρ P p ρ -1 α p . Hence P p ρ -1 α p -ρ -1 P p α p L 2 (Ω) = ρ -1 P p ρ P p ρ -1 α p -ρ -1 P p α p L 2 (Ω)

≤ 1 √ ρ - 1 √ ρ P p ρP p ρ -1 α p -α p L 2 (Ω) ≤ 1 √ ρ - ρP p ρ -1 α p -ρρ -1 α p L 2 (Ω) ≤ ρ + √ ρ - P p ρ -1 α p -ρ -1 α p L 2 (Ω) ,
where we have used the lower bound ρ -for ρ from the first to the second line. From the second to the third line, we have used the norm of the projector P p , which is 1 thanks to its orthogonality in L 2 (Ω) endowed with the inner product v 1 , ρ -1 v 2 L 2 , for all v 1 , v 2 ∈ L 2 (Ω). Finally, we have used the upper bound ρ + for ρ from the third to the fourth line. By (H1), still with v p = ρ -1 α p , we get

P p ρ -1 α p -ρ -1 P p α p L 2 (Ω) ≤ ρ + C p √ ρ - h θp ρ -1 α p H κ+1 (Ω) ,
leading to the announced result.

ISAE-SUPAERO, Université de Toulouse, France E-mail: ghislain.haine@isae.fr and denis.matignon@isae.fr and serhani@cerfacs.fr

Remark 2 . 3 .

 23 The compatibility condition u(0) = γ ⊥ T ---→ grad(w 0 ) is well-known in boundary control systems theory. See for instance [55, Chapter 10] for more details.

  thus estimates (H1)-(H2)-(H3) can indeed be applied to -→ α q and α p for all time t ∈ [0, T ].

Corollary 4 . 3 .

 43 Under the assumptions of Theorem 4.2, it holds

Figure 2 .

 2 Figure 2. Convergence of E X as a function of h at order 2, with the optimal choice of finite element families. See values in bold on Table2.

Figure 3 .

 3 Figure 3. Convergence of E X as a function of h at order 3, with the optimal choice of finite element families. See values in bold on Table3.

DG0

  boundary finite elements on ∂Ω

6. 4 . 1 .

 41 A non-convex case: the L-shaped domain. In this section, we test a non-convex case, with an analytical solution given as in the previous tests (boundary control is again given by the restriction of the analytical solution to ∂Ω), with a time step dt = 10 -3 . One can appreciate on Figure4how our results carry over, and remain valid even for this more complicated geometric shape.

Figure 4 .

 4 Figure 4. L-shape: a non-convex domain.

Figure 5 .

 5 Figure 5. An anisotropic case.

5 1. 5 ,

 55 For our example on Ω := -→ x ∈ R 2 | -→ x < 1 , the non-uniform anisotropic elasticity tensor T and the non-uniform heterogeneous mass density ρ are taken as followsT (x, y) := 2 0.2(1 + x)(1 -x) 0.2(1 + x)(1 -x) 1 , ρ(x, y) := 2 + 0.25(1 + x)(1 -x), ∀(x, y) ∈ Ω.The boundary control is chosen as u(t, x, y) := 5x sin(t) sin(1 -t), ∀t < 1, (x, y) ∈ ∂Ω, 0, ∀t > 1, (x, y) ∈ ∂Ω, while the admittance is defined by Y (t, x, y) := 2.5x sin(t) sin t-1.∀t > 1.5, (x, y) ∈ ∂Ω, 0, ∀t < 1.5, (x, y) ∈ ∂Ω,

Figure 6 .

 6 Figure 6. Time evolution of the energies present in the system (RT 1 × CG 1 × CG 1 discretization). The number of degrees of freedom is: 2055 + 714 + 84 = 2853.

L 2

 2 grad ρ -1 α p -P p ρ -1 α p L 2 (Ω) + ---→ grad P p ρ -1 α p -ρ -1 P p α p L 2 (Ω)

  The proposed strategy can handle other causalities, i.e. other collocated boundary control and observation, in a straightforward manner.

		e q (t)	
	T  on the left leads to -1 d d t e p (t)  u(t) d d t e q (t) + e p (t) M ρ d d t e p (t) = (u(t)) M ∂ y(t), thanks to the skew-symmetry of the extended structure matrix, and the result e q (t) M T -1 follows. e Multiplying (15) by 3.2. Other causalities.

is straightforward. Clearly d dt H d (t) = e q (t) M q (t) + e p (t) M ρ d d t e p (t) , thanks to the symmetry of the mass matrices. Let us focus on the other uniform causality, i.e. with deflection velocity control. It only consists on switching the role played by u and y, i.e. replace (2) by (2S)

Table 2 .

 2 Convergence at order 2 for E X

Table 4 .

 4 Hamiltonian convergence rates for DG0 boundary finite elements. The possible optimal order given by Theorem 4.4 is θ * = 2 (i.e. twice the order reached for the state error convergence rate shown in Table1) cannot be reached.

Table 5 .

 5 Hamiltonian convergence rates for DG1 and CG1 boundary finite elements. The possible optimal order given by Theorem 4.4 is θ * = 2 or 4, i.e. twice the orders reached for the state error convergence rate shown in Table2. In bold the convergence rate achieving order 2 while state error is only of order 1, i.e. a numerical evidence that compatibility conditions must be satisfied in those cases.

	0.82 0.81 0.8	0.83 0.82 0.81	0.84 0.81 0.81
	0.8	0.81	0.81
	0.8 0.8 0.81 0.8 0.81 0.81	0.81 0.66 0.82 0.81 0.81 0.82	0.81 1.09 0.99 0.81 0.82 0.81	DG0 boundary finite elements on ∂Ω
	0.81	0.82	0.81
	8	0.81	0.8
		-0.1	-0.1
		CG2	CG3

Table 6 .

 6 Hamiltonian convergence rates for DG2 and CG2 boundary finite elements. The possible optimal order given by Theorem 4.4 is θ * = 2, 4 or 6, i.e. twice the orders reached for the state error convergence rate shown in Table3. In bold the convergence rate achieving order 2 (resp. 4) while state error is only of order 1 (resp. 2), i.e. a numerical evidence that compatibility conditions must be satisfied in those cases.

	2.09 2.08 2.07	2.04 2.02 2.02	1.14 2.02 2.03
	2.07	2.02	2.02
	2.07 1.96 2.08 2.07 2.04 2.08	2.02 1.96 2.02 2.02 2.06 2.02	2.02 2.11 2.02 2.03 2.0 2.02	CG1 boundary finite elements on ∂Ω
	2.07	2.02	2.02
	1.99	1.95	1.99
	2.05	0.01	0.04
		CG2	CG3

  ,

	the kinetic energy
	E Kin (t) :=	1 2 Ω	α q (t, -→ x ) ρ( -→ x )	2	d -→ x ,
	the boundary-supplied energy
	S(t) :=	0	t	u(t, -→ s ), y(t, -→ s ) U ,Y dt,
	and the damped energy
	D(t) :=	0	t	Y (t, -→ s ) y(t, -→ s ), y(t, -→ s ) U ,Y dt.

The total energy present in the system over time is then decomposed as

E(t) = E Pot (t) + E Kin (t) + S(t) + D(t) = H(t) + S(t) + D(t),

  ∂ t P p α p -∂ t α d p , ρ -1 P p α p -α d p L 2 (Ω) ,leading to the announced result.Lemma A.3. Under the assumptions of Theorem 4.2, one has for all h ∈ (0, h

* ) ---→ grad ρ -1 (α p -P p α p ) L 2 (Ω)

Acknowledgments

This work has been performed in the frame of the Collaborative Research DFG and ANR project INFIDHEM, entitled Interconnected Infinite-Dimensional systems for Heterogeneous Media, n • ANR-16-CE92-0028. Further information is available at https://websites.isae-supaero.fr/infidhem/the-project/

Appendix A. Technical lemmas

The following lemma gives an upper bound for θ 1,0 in (H5), namely θ 0,1 from (H4), although not optimal for general geometry (e.g. for convex domains).

Lemma A.1. If (H4) holds true, then there exists a constant C 0,1 > 0 such that

where we have used P p P 1,p v p = P 1,p v p for all v p ∈ H 1 (Ω), using (H4), P p L(L 2 (Ω)) = 1, and we have defined

We provide here the two technical lemmas used in the proof of Theorem 4.2.

Lemma A.2. Under the assumptions of Theorem 4.2, one has

Proof. From the weak formulations ( 10)- [START_REF] Brugnoli | Partitioned Finite Element Method for Power-Preserving Structured Discretization with Mixed Boundary Conditions[END_REF], and thanks to the conformity of the finite element families, one has

and

Summing the latter two equalities gives

---→ grad ρ -1 P p α p -α d p L 2 (Ω)

+ u -u d , γ 0 ρ -1 P p α p -α d p U ,Y .