P Kotyczka

B Maschke

L Lefèvre

Weak Form of the Stokes-Dirac Structure and Geometric Discretization of Port-Hamiltonian Systems

Keywords: Open systems of conservation laws, port-Hamiltonian systems, mixed Galerkin methods, geometric discretization, structure-preserving discretization

We present the mixed Galerkin discretization of distributed-parameter port-Hamiltonian systems. Due to the inherent definition of (boundary) interconnection ports, this system representation is particularly useful for the modeling of interconnected multi-physics systems and control. At the prototypical example of a system of two conservation laws in arbitrary spatial dimension, we derive the main contributions: (i) A weak formulation of the underlying geometric (Stokes-Dirac) structure, (ii) its geometric approximation by a finite-dimensional Dirac structure using a mixed Galerkin approach and power-preserving maps on the space of discrete power variables and (iii) the approximation of the Hamiltonian to obtain finite-dimensional port-Hamiltonian state space models. The powerpreserving maps on the discrete bond space offer design degrees of freedom for the discretization, which is illustrated at the example Whitney finite elements on a 2D simplicial triangulation. The resulting schemes can be considered as trade-offs between centered approximations and upwinding.

Introduction

The port-Hamiltonian (PH) approach for the modeling, interconnection and control of multi-physics systems underwent an enormous evolution during the past two decades. In this article, we concentrate on distributed-parameter PH systems as initially presented in [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF], and refer the reader to the books [START_REF] Duindam | Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach[END_REF], [START_REF] Jacob | Linear port-Hamiltonian systems on infinite-dimensional spaces[END_REF] and [START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF] for a more general overview on theory and applications. The salient feature of PH systems is their representation in terms of (i) a linear geometric structure -a Stokes-Dirac structure -that describes the power flows inside the system and over its boundary and (ii) an energy functional (or more generally a potential) from which the constitutive or closure relations are derived, and which determines the nature of the systems. Hence, completely different systems -linear/nonlinear or hyperbolic/parabolic [START_REF] Zwart | Building systems from simple hyperbolic ones[END_REF] can share the same interconnection structure. PH systems are by definition open systems, they interact with their environment through boundary ports. The in-and outputs in the sense of systems' theory and control are defined via a duality product whose value equals the exchanged power at the port. The appropriate definition of boundary port variables plays a crucial role in showing that a PH system is a well-posed boundary control systems [START_REF] Le Gorrec | Dirac structures and boundary control systems associated with skew-symmetric differential operators[END_REF]. The reasonable definition of distributed power variables as inand outputs is discussed in [START_REF] Nishida | Boundary integrability of multiple stokes-dirac structures[END_REF].

The simulation and control by numerical methods, of complex (in the sense of their geometry, their nonlinearities, their extent, their couplings) distributed-parameter PH systems, requires a spatial discretization, which shall retain the underlying geometric properties related to power continuity. According to the separation of the interconnection structure from the constitutive equations, a geometric discretization consists of two steps:

• Finite-dimensional approximation of the underlying Stokes-Dirac structure. The duality between the power variables variables (their duality product has the interpretation of power) must be mapped onto the finite-dimensional approximation. This requires a mixed approach with different approximation spaces for each group of dual power variables (called flows and efforts). The subspace of discrete power variables on which (the preserved) power-continuity holds, defines a Dirac structure as a finite-dimensional counterpart of the Stokes-Dirac structure.

• The discretization of the interconnection structure is followed by an appropriate discretization of the Hamiltonian and the closure equations, again taking into account the different nature of the approximation spaces.

A geometric discretization is, hence, a compatible discretization as defined in [START_REF] Pavel | Principles of mimetic discretizations of differential operators[END_REF]: "Compatible discretizations transform partial differential equations to discrete algebraic problems that mimic fundamental properties of the continuum equations." The open character of PH systems requires special attention to the treatment of the boundary port variables, in particular the boundary inputs which are imposed as boundary conditions.

The first approach for a geometric or structure-preserving discretization of PH systems [START_REF] Golo | Hamiltonian discretization of boundary control systems[END_REF], application to a diffusive process in [START_REF] Baaiu | Structure-preserving infinite dimensional model reduction: Application to adsorption processes[END_REF], is based on the idea of mixed finite elements. The Stokes-Dirac structure is, however, discretized in strong form which produces restrictive compatibility conditions. In the 1D pseudo-spectral method [START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF], the degeneracy of the discrete duality product is rectified by the definition of reduced effort variables, see also the recent paper [START_REF] Minh | Symplectic spatial integration schemes for systems of balance equations[END_REF] for the application to plasma dynamics described by a parabolic PDE. The reduction of the discrete effort space is generalized by the power-preserving maps in the present paper. In [? ], the Stokes-Dirac structure is reformulated, changing the role of state and co-state variable for one conservation law. The discrete power-variables are immediately connected with a non-degenerate duality pairing, at the price of a metric-dependent interconnection structure. The simplicial discretization of PH systems based on discrete exterior calculus has been introduced in [START_REF] Seslija | Discrete exterior geometry approach to structure-preserving discretization of distributed-parameter port-hamiltonian systems[END_REF], see [START_REF] Kotyczka | Discrete port-Hamiltonian formulation and numerical approximation for systems of two conservation laws[END_REF] for the extended definition of boundary inputs and a revised interpretation of the resulting state space models. Recently, the very related discretization on staggered grids has been reported using finite volumes [START_REF] Kotyczka | Finite volume structure-preserving discretization of 1D distributed-parameter port-Hamiltonian systems[END_REF] and finite differences [START_REF] Trenchant | Structure preserving spatial discretization of 2D hyperbolic systems using staggered grids finite difference[END_REF].

The weak or variational formulation as the basis for Galerkin numerical approximations, including the different varieties of the finite element method (see [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF], to cite only one textbook), has been only rarely used for modeling and discretization of PH systems: In [? ], one of the two conservation laws is written in weak form. [START_REF] Altmann | A port-hamiltonian formulation of the navier-stokes equations for reactive flows[END_REF] presents the PH model of the reactive 1D Navier-Stokes equations in weak form. In [START_REF] Luiz Cardoso-Ribeiro | A port-Hamiltonian model of liquid sloshing in moving containers and application to a fluid-structure system[END_REF], a piezo patch on a flexible beam, with a discontinuous spatial distribution, is included in the PH model and its structure preserving discretization via the weak form.

In this article, we present the geometric discretization of distributed-parameter PH systems based on the weak or variational formulation of the underlying Stokes-Dirac structure. Doing so, some limitations and restrictions of current approaches for PH systems can be overcome.

• The formulation is valid for systems on spatial domains with arbitrary dimension.

• PH systems with non-canonical system operators (containing e.g. higher order derivatives) can be discretized, as well as systems with dissipation and/or diffusion.

• In terms of finite elements, approximation spaces with higher polynomial degree can be used.

• Due to the weak imposition of boundary conditions, the boundary inputs appear directly in the weak form of the Stokes-Dirac structure.

• The power-preserving maps for the discrete power variables offer design degrees of freedom for the resulting PH approximate models. These design parameters can be used to adapt the discretization to the physical nature of the system (e.g. accounting for the ratio between convection and diffusion) and to improve the approximation quality.

We consider as the prototypical example of distributed-parameter PH systems, an open system of two conservation laws in canonical form, as presented in [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF]. We use the language of differential forms, see e.g. [START_REF] Flanders | Differential Forms with Applications to the Physical Sciences[END_REF], which highlights the geometric nature of each variable and allows for a unifying representation independent of the dimension of the spatial domain.

An important reason for expressing the spatial discretization of PH systems based on the weak form is to make the link with modern geometric discretization methods. Bossavit's work in computational electromagnetism, [START_REF] Bossavit | Differential forms and the computation of fields and forces in electromagnetism[END_REF], [START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements[END_REF] (based on Whitney elements [START_REF] Whitney | Geometric integration theory[END_REF]) and Tonti's cell method [START_REF] Tonti | A direct discrete formulation of field laws: The cell method[END_REF] keep track of the geometric nature of the system variables which allows for a direct interpretation of the discrete variables in terms of integral system quantities. This integral point of view is also adopted in discrete exterior calculus [START_REF] Mathieu Desbrun | Discrete differential forms for computational modeling[END_REF]. Finite element exterior calculus [START_REF] Arnold | Finite element exterior calculus: from hodge theory to numerical stability[END_REF] gives a theoretical frame to describe functional spaces of differential forms and their compatible approximations, which includes the construction of higher order approximation bases that generalize the famous Whitney forms [START_REF] Whitney | Geometric integration theory[END_REF], see also [START_REF] Rapetti | Whitney forms of higher degree[END_REF]. We refer also to the recent article [START_REF] Rr Hiemstra | High order geometric methods with exact conservation properties[END_REF] which proposes conforming polynomial approximation bases, in which the conservation laws are exactly satisfied, and which gives an excellent introduction to the geometric discretization. Impressing examples for the use of geometric discretization methods can be found in weather prediction [START_REF] Cotter | Mixed finite elements for numerical weather prediction[END_REF] or the simulation of large-scale fluid flows [START_REF] Cotter | A finite element exterior calculus framework for the rotating shallow-water equations[END_REF], where the conservation of potential vorticity plays an important role. Another important aspect of using the weak form as basis for structurepreserving discretization is to make the link with well-known numerical methods and to pave the way for a simulation of PH systems with existing numerical tools like FreeFEM++ [START_REF] Hecht | New development in freefem++[END_REF] or GetDP [START_REF] Geuzaine | GetDP: a general finite-element solver for the de Rham complex[END_REF].

The paper is structured as follows. In Section 2 we review the definition of distributed-parameter port-Hamiltonian systems. We focus in particular on systems of two conservation laws and the underlying Stokes-Dirac structure. We introduce the Stokes-Dirac structure with boundary ports of different causality as it will be used in the rest of the paper. In Section 3 we define the weak form of this Stokes-Dirac structure and present its mixed Galerkin approximation. This approximation features, due to the different geometric nature of the power variables and their approximation spaces, a degenerate power pairing.

In Section 4, new discrete power variables (or bond variables), with a non-degenerate duality pairing are defined via power-preserving maps. The mixed Galerkin approximation, expressed in the new space of discrete bond variables, defines a finite-dimensional Dirac structure, which admits different representations. The input-output representation, together with the finite-dimensional approximation of the Hamiltonian, leads to the desired port-Hamiltonian approximate models in state space form. Section 5 illustrates the method using Whitney finite elements on a 2D simplicial grid. In particular, we highlight the interpretation of the finite-dimensional state and power variables in terms of integral quantities on the grid. Moreover, we show how discretization schemes with different approximation quality are obtained via the parametrization of the power-preserving maps. Section 6 closes the paper with a summary and an outlook to ongoing and future work.

2 Distributed-Parameter Port-Hamiltonian Systems

Canonical systems of two conservation laws

We consider canonical systems of two conservation laws 1 on an open, bounded and connected n-dimensional spatial domain 2 Ω with Lipschitz boundary ∂Ω in port-Hamiltonian form as described in [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF]. We use the language of differential forms, see e.g. [START_REF] Flanders | Differential Forms with Applications to the Physical Sciences[END_REF]. Initially, we assume all differential forms to be sufficiently smooth and denote Λ k (Ω) the space of smooth differential k-forms on Ω. The partial differential equations (PDEs) for this canonical class of PH systems are written as 3

-∂ t p(z, t) -∂ t q(z, t) = 0 (-1) r d d 0 Jcan δ p H(p(z, t), q(z, t)) δ q H(p(z, t), q(z, t)) , (1)
with p(z, t) ∈ Λ p (Ω), q(z, t) ∈ Λ q (Ω) the state variables subject to initial conditions p(z, 0) = p 0 (z), q(z, 0) = q 0 (z). J can denotes a canonical, matrix-valued and formally skew-symmetric 4 differential operator with d : Λ k (Ω) → Λ k+1 (Ω) the exterior derivative. Moreover, r = pq + 1 5 . The Hamiltonian or energy functional is given by

H(p(z, t), q(z, t)) = Ω H(p(z, t), q(z, t), z) (2)
with the Hamiltonian density H : Λ p (Ω) × Λ q (Ω) × Ω → Λ n (Ω). The vector on the right hand side of (1) contains the co-energy or co-state variables e p (z, t) := δ p H(p(z, t), q(z, t)) ∈ Λ n-p (Ω), e q (z, t) := δ q H(p(z, t), q(z, t)) ∈ Λ n-q (Ω),

which are -for the considered class of hyperbolic systems -defined as the variational derivatives of the Hamiltonian, i.e. differential forms δ p H and δ q H that satisfy6

H(p + δp, q + δq) = Ω H(p, q, z) + Ω δ p H ∧ δp + δ q H ∧ δq + O(δp 2 , δq 2 ). ( 4
)
We call the equations to define e p and e q constitutive or closure equations. They determine the nature of the system (linear or nonlinear; hyperbolic as considered here or parabolic, see [START_REF] Zwart | Building systems from simple hyperbolic ones[END_REF]). Comparing the degrees of the co-energy differential forms (np and nq) with the degrees for the right hand side of (1) to make sense (q -1 and p -1) yields the relation p + q = n + 1 between the dimension of Ω and the degrees of the differential forms 7 .

The boundary port variables are defined as

f ∂ (z, t) e ∂ (z, t) = 1 0 0 (-1) p δ p H(p(z, t), q(z, t)) δ q H(p(z, t), q(z, t)) ∂Ω , (5)
where (•)| ∂Ω represents the trace or the restriction to the boundary of a differential form. Imposing the port variables f ∂ (z, t) ∈ Λ n-p (∂Ω) and/or e ∂ (z, t) ∈ Λ n-q (∂Ω) on a subset of ∂Ω as control input (and understanding the port variables on the remaining boundary as observation or output), makes the system representation (1), ( 5) a boundary control system in the sense of [START_REF] Hector | Boundary control systems[END_REF].

For 1D linear PH systems with a generalized skew-symmetric system operator, [START_REF] Le Gorrec | Dirac structures and boundary control systems associated with skew-symmetric differential operators[END_REF] gives conditions on the assignment of boundary in-and outputs for the system operator to generate a contraction semigroup. The latter is instrumental to show well-posedness of a linear PH system, see [START_REF] Jacob | Linear port-Hamiltonian systems on infinite-dimensional spaces[END_REF]. Essentially, at most half the number of boundary port variables can be imposed as control inputs for a well-posed PH system in 1D.

Taking δp = ṗ, δq = q as variations in (4), and omitting the higher order terms, the time derivative of the energy functional (2) reads8

Ḣ = Ω δ p H ∧ ṗ + δ q H ∧ q =: δ p H| ṗ Ω + δ q H| q Ω , (6)
Replacing ṗ, q according to (1), using the integration-by-parts formula for differential forms and the definition (5) of boundary port variables, yields

Ḣ = ∂Ω (-1) p δ q H| ∂Ω ∧ δ p H| ∂Ω =: (-1) p δ q H|δ p H ∂Ω . ( 7
)
Equating the right hand sides of the last two equations gives a (differential) energy balance equation or a power continuity equation that relates the internal power flows (i.e. the lossless exchange of energy between the two energy types) inside the domain Ω (6) and the power supplied over the system boundary ∂Ω [START_REF] Baaiu | Structure-preserving infinite dimensional model reduction: Application to adsorption processes[END_REF]. Due to the energy exchange over the system boundary, we call (1), [START_REF] Arnold | Finite element exterior calculus: from hodge theory to numerical stability[END_REF] an open system of two conservation laws in canonical port-Hamiltonian form.

Notation 1. •|• Ω : Λ n-k (Ω) × Λ k (Ω) → R
is the short notation for the duality product

λ|µ Ω := Ω λ ∧ µ (8)
between two differential forms on the n-dimensional domain Ω. Accordingly for ∂Ω.

Example 1 (1D transmission line). The simplest example in 1D of a system of two conservation laws is an electric transmission line with Ω = [0, L], see e.g. [START_REF] Golo | Approximation of the telegrapher's equations[END_REF]. With p(z) ∈ Λ 1 (Ω), the magnetic flux density one-form, q(z) ∈ Λ 1 (Ω), the electric charge density one-form, l(z)dz, c(z)dz ∈ Λ 1 (Ω) the distributed inductance and capacitance per length (note that in this representation, l(z) and c(z) are functions, i.e. zero-forms), the Hamiltonian density one-form is H(p, q) = 1 2 p(z) ∧ * p(z) l(z) + q(z) ∧ * q(z) c(z) . The Hodge star operator * : Λ k (Ω) → Λ n-k (Ω) renders in the 1D case a one-form a zero-form and vice-versa 9 . The variational derivatives of the Hamiltonian H = L 0 H turn out to be the current and the voltage along the line,

δ p H = * p(z) l(z) = i(z) ∈ Λ 0 (Ω), δ q H = * q(z) c(z) = u(z) ∈ Λ 0 (Ω).
Example 2 (2D shallow water equations). The shallow water equations describe the two-dimensional flow of an inviscid fluid with relatively low depth ("shallow") which permits the averaging of the horizontal components of the velocity field and the omission of the vertical velocity component. The two equations that describe the conservation of mass and momentum over an infinitesimal, fixed surface element 10 , are (we consider the fluid in a non-rotating system) can be written in vector calculus notation, see e.g. [START_REF] Fernandez-Nieto | Coupling superposed 1d and 2d shallow-water models: Source terms in finite volume schemes[END_REF],

∂ t h + div(hu) = 0 ∂ t (hu) + div(huu) + 1 2 g∇h 2 + gh∇z b = 0. ( 10
)
h denotes the water level over the bottom, z b is the vertical coordinate that describes the bottom profile, u = [u v] T the 2-dimensional velocity field, hu = F the discharge vector and g the gravitational acceleration. u • u and uu denote the scalar and the tensor (dyadic) product of two vectors, respectively. With some rules of tensor calculus 11 , and replacing the continuity equation, the momentum equation can be reformulated in terms of u and one obtains

∂ t h ∂ t u + qF ⊥ = 0 -div -grad 0 1 2 u • u + gh + gz b hu , ( 11
)
where q = 1 h (∂ x v -∂ y u) denotes the potential vorticity 12 , and F ⊥ = [hvhu] T . The term qF ⊥ represents the acceleration of the fluid due to rotation of the flow. It stems from the rotational part of the transport term in the momentum equation. The total energy (per unit mass) is

H = Ω 1 2 hu • u + 1 2 gh 2 + ghz b dx. ( 12
)
To rewrite the equations in terms of differential forms (in a covariant formulation), we use the relations, see [START_REF] Abraham | Manifolds, tensor analysis, and applications[END_REF]

13 ∇f = (df ) , div f = * d( * f ). ( 13
)
9

The Hodge star induces an inner product on the space of differential forms on a manifold M by

(α, β) := M α ∧ * β = M β ∧ * α = (β, α), α, β ∈ Λ k (M ). (9)
see Section 8.4 of [START_REF] Flanders | Differential Forms with Applications to the Physical Sciences[END_REF] or Section 3.6 of [START_REF] Holm | Geometric mechanics: Dynamics and symmetry[END_REF]. The inner product is not necessarily the standard L 2 norm, but may be equipped with another metric, see e.g. the energy norm fo linear PH systems [START_REF] Jacob | Linear port-Hamiltonian systems on infinite-dimensional spaces[END_REF]. The Hodge star is, hence, metric dependent. 10 Which corresponds to the Eulerian representation of the fluid flow.

11 See appendix A.4 of [START_REF] Bird | Transport phenomena[END_REF]:

∇ • (sI) = ∇s, ∇ • (vw) = v • ∇w + w(∇ • v), v • ∇v = 1 2 ∇(v • v) -v × (∇ × v).
The last term with cross product and rotation has to be evaluated based on the 3D velocity vector with zero vertical component. 12 The potential vorticity satisfies the balance equation ∂tq + u • ∇q = 0, i.e. it is advected with the fluid flow see e.g. [START_REF] Arakawa | A potential enstrophy and energy conserving scheme for the shallow water equations[END_REF]. It plays an important role in the long-time numerical simulation of large scale flow problems, see e.g. [START_REF] Td Ringler | A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured c-grids[END_REF]. 13 Index raising ( ) produces a vector field with the same components from a one-form. Index lowering ( ) produces a one-form with identical components from a vector field. Raising and lowering refers to the fact that upper (lower) indices are typically used for the components of vector fields (one-forms).

Taking into account that * * λ = (-1) k(n-k) λ for a k-form λ, we obtain

-∂ t ( * h) -(∂ t u + qF ⊥ ) = 0 -d d 0 p dyn -( * F ) , ( 14
)
where * h ∈ Λ 2 (Ω) and u ∈ Λ 1 (Ω) are the 2-form and 1-form associated with the water depth and the flow velocity (p = 2, q = 1). p dyn = 1 2 u • u + gh + gz b ∈ Λ 0 (Ω) is the hydrodynamic pressure function (0-form) and * F ∈ Λ 1 (Ω) is the 1-form associated to the discharge per unit width. Indeed the vector on the right can be expressed in terms of the variational derivatives

p dyn = δ * h H and -( * F ) = δ u H of the Hamiltonian H = Ω H density 2-form 14 H = 1 2 h(u ∧ * u ) + 1 2 gh * h + g * hz b .
If the rotational term qF ⊥ can be neglected15 , (14) has the canonical structure (1).

Stokes-Dirac structures and port-Hamiltonian systems

Power continuity is a property of the underlying linear geometric structure that relates the dual power variables of the system, independent of the closure equations for the co-state variables. Defining the distributed and boundary flow and effort differential forms as

f p f q := -∂ t p -∂ t q , e p e q := δ p H δ q H , f ∂ e ∂ := δ p H (-1) p δ q H ∂Ω , ( 15
)
the Stokes-Dirac structure describes a linear subspace of the bond space 16 F × E,

F = Λ p (Ω) × Λ q (Ω) × Λ n-p (∂Ω) E = Λ n-p (Ω) × Λ n-q (Ω) × Λ n-q (∂Ω) (16)
of dual power variables where power continuity e p |f p Ω + e q |f q Ω power extracted from distributed storage + e ∂ |f ∂ ∂Ω power supplied over the boundary = 0 [START_REF] Luiz Cardoso-Ribeiro | A port-Hamiltonian model of liquid sloshing in moving containers and application to a fluid-structure system[END_REF] holds. The Stokes-Dirac structure extends the notion of the finite-dimensional Dirac structure to infinitedimensional systems. The latter can be considered "the geometrical notion formalizing general powerconserving interconnections" [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF]:

Definition 1 ([60], Def. 2.1
). Given the finite-dimensional linear spaces F fin and E fin over R17 and a duality pairing

•|• : F fin × E fin → R. Define the symmetric bilinear form (f 1 , e 1 ), (f 2 , e 2 ) := e 1 |f 2 + e 2 |f 1 , (f i , e i ) ∈ F × E, i = 1, 2. (18)
A Dirac structure is a linear subspace D fin ⊂ F fin × E fin such that D ⊥ fin = D fin , where ⊥ denotes the orthogonal complement with respect to the bilinear form •, • .

From the requirement that each element of a Dirac structure D fin belongs to its orthogonal complement

D ⊥
fin (and vice versa), it follows that for all (f , e) ∈ D fin

(f , e), (f , e) = 2 e|f = 0, (19)
which is a (finite-dimensional) power-continuity relation as [START_REF] Luiz Cardoso-Ribeiro | A port-Hamiltonian model of liquid sloshing in moving containers and application to a fluid-structure system[END_REF]. For more details and different representations of Dirac structures (in the PH context) refer to [START_REF] Van Der Schaft | L2-Gain and Passivity Techniques in Nonlinear Control[END_REF], [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF].

For the infinite-dimensional Stokes-Dirac structure, we consider the duality pairing over differential forms as defined in [START_REF] Bird | Transport phenomena[END_REF] and the symmetric bilinear form

(f p 1 , f q 1 , f ∂ 1 , e p 1 , e q 1 , e ∂ 1 ), (f p 2 , f q 2 , f ∂ 2 , e p 2 , e q 2 , e ∂ 2 ) := e p 1 |f p 2 Ω + e q 1 |f q 2 Ω + e ∂ 1 |f ∂ 2 ∂Ω + e p 2 |f p 1 Ω + e q 2 |f q 1 Ω + e ∂ 2 |f ∂ 1 ∂Ω (20)
with (i = 1, 2)

f p i ∈ Λ p (Ω), f q i ∈ Λ q (Ω), f ∂ i ∈ Λ n-p (Ω) e p i ∈ Λ n-p (Ω), e q i ∈ Λ n-q (Ω), e ∂ i ∈ Λ n-q (Ω). ( 21
)
Theorem 1. ( [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF], Thm. 2.1) The linear subspace of F × E, see [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF], given by

D = {(f p , f q , f ∂ , e p , e q , e ∂ ) ∈ F × E | f p f q = 0 (-1) r d d 0 e p e q , f ∂ e ∂ = 1 0 0 (-1) p e p | ∂Ω e q | ∂Ω }. ( 22
)
with r = pq + 1 and | ∂Ω the restriction to the boundary ∂Ω is a Dirac structure, i.e. D = D ⊥ holds with respect to the bilinear mapping defined in [START_REF] Cotter | A finite element exterior calculus framework for the rotating shallow-water equations[END_REF].

As the application of Stokes' theorem instrumental in showing that the above subspace is a Dirac structure, the so-defined geometric structure is called a Stokes-Dirac structure. In order to represent the situation of boundary inputs with different causality in a formally consistent way, we introduce the Stokes-Dirac structure with an alternative designation of boundary port variables.

Corollary 1. The linear subspace

D = {(f p , f q , ê∂ , e p , e q , f ∂ ) ∈ F × Ê | f p f q = 0 (-1) r d d 0 e p e q , ê∂ f ∂ = 1 0 0 (-1) p e p | ∂Ω e q | ∂Ω } (23) with F = Λ p (Ω) × Λ q (Ω) × Λ n-q (∂Ω) Ê = Λ n-p (Ω) × Λ n-q (Ω) × Λ n-p (∂Ω) ( 24
)
is a Dirac structure.

We can understand the representation of the Hamiltonian systems with canonical interdomain coupling (1), ( 5) as the composition of an underlying geometric interconnection for the pairs of port variables (f p , e p ), (f q , e q ), (f ∂ , e ∂ ) with dynamic and closure equations [START_REF] Brenner | The mathematical theory of finite element methods[END_REF].

Definition 2. ([60], Def. 2.2, Prop. 2.1)
The distributed-parameter system on a closed and bounded ndimensional manifold Ω with boundary ∂Ω, defined by ( 1), [START_REF] Arnold | Finite element exterior calculus: from hodge theory to numerical stability[END_REF] with Hamiltonian H(p, q) = Ω H(p, q, z) and Hamiltonian density

H : Λ p (Ω) × Λ q (Ω) × Ω → Λ n (Ω), with p + q = n + 1 and r = pq + 1 is a port-Hamiltonian system.
It represents an open system of two conservation laws and satisfies the energy balance equation [START_REF] Baaiu | Structure-preserving infinite dimensional model reduction: Application to adsorption processes[END_REF].

Boundary ports

The choice of boundary port variables to define a Stokes-Dirac structure (an infinite-dimensional PH system) is not unique, see [START_REF] Le Gorrec | Dirac structures and boundary control systems associated with skew-symmetric differential operators[END_REF] for the 1D case. On parts of the boundary, e q may define the control input, while, vice versa, this role is assigned to e p on the rest. The only constraint on the definition of pairs of boundary port variables is that their entirety reflects the power flow over the boundary as indicated in [START_REF] Baaiu | Structure-preserving infinite dimensional model reduction: Application to adsorption processes[END_REF].

In order to represent a larger class of boundary control problems for systems of two conservation laws, the following proposition generalizes the definition of the Stokes-Dirac structure to multiple pairs of inand outputs on ∂Ω with different causality.

f Γ i = e p | Γi , e Γ i = (-1) p e q | Γi , f Γ j = (-1) p e q | Γj , êΓ j = e p | Γj , ( 25
)
such that the bond space

F × E is composed of 18 F = Λ p Ω × Λ q Ω × Λ n-p Γ1 × • • • × Λ n-p Γn Γ × Λ n-q Γ1 × • • • × Λ n-q Γn Γ E = Λ n-p Ω × Λ n-q Ω × Λ n-q Γ1 × • • • × Λ n-q Γn Γ × Λ n-p Γ1 × • • • × Λ n-p Γn Γ . ( 26
)
The subspace D ⊂ F × E, on which

f p f q = 0 (-1) r d d 0 e p e q ( 27
)
holds and the boundary ports are defined as in [START_REF] Fernandez-Nieto | Coupling superposed 1d and 2d shallow-water models: Source terms in finite volume schemes[END_REF], is a Dirac structure.

Proof. First observe that with the choice of boundary ports, and by construction of the subsets Γ i and Γj , the boundary power flow can be expressed as

n Γ i=1 e Γ i |f Γ i Γi + nΓ j=1 f Γ j |ê Γ j Γj = (-1) p e q |e p ∂Ω . ( 28
)
The proof that the above subspace is a Dirac structure, consists of appropriately decomposing Ω and exploiting the compositionality property, see Remark 2.2 of [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF], of the Stokes-Dirac structure on each subset.

1. Decompose Ω in a set of n-dimensional submanifolds Ω k and Ωl , with the same orientation as ∂Ω on ∂Ω k ∩ ∂Ω and ∂ Ωl ∩ ∂Ω. On each subset, a Stokes-Dirac structure is defined, in an alternating manner, according to Theorem 1 and Corollary 1, respectively. Then Remark 2. In the above proposition, boundary efforts and flows are defined as pure restrictions of either of the distributed efforts to the corresponding subsets Γ i , Γj of the boundary. It is, however, also possible to define mixed boundary variables (e.g. scattering variables), as long as the boundary power continuity equation [START_REF] Geuzaine | GetDP: a general finite-element solver for the de Rham complex[END_REF] holds. Such a more general parametrization of boundary in-and outputs is presented in [START_REF] Le Gorrec | Dirac structures and boundary control systems associated with skew-symmetric differential operators[END_REF] for linear 1D PH systems with general skew-symmetric interconnection operators. In-and outputs are characterized such that the resulting system is a boundary control system which is associated to a contraction semigroup.

k ∂Ω k ∩ ∂Ω = i Γ i , l ∂ Ωl ∩ ∂Ω = j Γj and Γ i ∩ ∂Ω kl = ∅, Γj ∩ ∂ Ωlk = ∅ for all i, j, k, l. ∂Ω kl = -∂
f p ∈ C 0 Λ p (Ω), f q ∈ C 0 Λ q (Ω), e p ∈ C 1 Λ n-p (Ω), e q ∈ C 1 Λ n-q (Ω). ( 29
)
The efforts being continuously differentiable on Ω, this holds also for their restriction to the boundary

∂Ω ⊂ Ω: e Γ i , f Γ j ∈ C 1 Λ n-q (∂Ω), êΓ j , f Γ i ∈ C 1 Λ n-p (∂Ω), i = 1, . . . , n Γ , j = 1, . . . , nΓ . ( 30
)
Note that with these functional spaces, the duality products (power density integrals) in ( 17) and ( 28) are well-defined and finite.

Further classes of distributed parameter PH systems

In this chapter, we concentrate on canonical systems of two conservation laws in arbitrary spatial dimension. Beyond this basic class of PH systems (which however covers different linear and nonlinear physical phenomena), there exists a growing number of PH models for different physical phenomena, see e.g. [START_REF] Vu | Port-hamiltonian formulation for systems of conservation laws: Application to plasma dynamics in tokamak reactors[END_REF] for the modeling of the plasma in a fusion reactor, [START_REF] Altmann | On the port-hamiltonian structure of the navier-stokes equations for reactive flows[END_REF] for the reactive Navier-Stokes flow or [START_REF] Ramirez | Irreversible port-hamiltonian systems: A general formulation of irreversible processes with application to the cstr[END_REF] for irreversible thermodynamic systems to mention only a few interesting examples. In [START_REF] Polner | A hamiltonian vorticity-dilatation formulation of the compressible euler equations[END_REF], a PH formulation of the compressible Euler equations in terms of density, weighted vorticity and dilatation is presented.

For the class of linear 1D PH systems represented by general higher order skew-symmetric system operators, the parametrization of in-and outputs to define a boundary control system with a system operator that generates a (contractive or unitary) C 0 semigroup is presented in [START_REF] Le Gorrec | Dirac structures and boundary control systems associated with skew-symmetric differential operators[END_REF]. The PH representation is not unique. An important approach for mechanical systems is based on a jet bundle formulation [START_REF] Schöberl | Jet bundle formulation of infinite-dimensional port-hamiltonian systems using differential operators[END_REF], see also [START_REF] Schoberl | Analysis and comparison of port-hamiltonian formulations for field theories-demonstrated by means of the mindlin plate[END_REF] where the PH formulation of the Mindlin plate based on the jet bundle approach and the Stokes Dirac structure are compared. In [START_REF] Schöberl | Differentialgeometrische beschreibung und analyse tor-basierter hamilton'scher systeme[END_REF], different PH representations (symplectic and polysymplectic, both based on jet bundles; Stokes-Dirac structure) which yield different definitions of state variables and Hamiltonian densities, are compared at the illustrative example of the Timoshenko beam.

Galerkin approximation of the Stokes-Dirac structure

In this section, we study the Galerkin approximation of the Stokes-Dirac structure defined in Prop. 1. To this end, we introduce the weak form of the Stokes-Dirac structure and use a mixed Galerkin approach to approximate the differential forms according to their geometric nature.

Weak form of the Stokes-Dirac structure

The first motivation to study the approximation of distributed-parameter PH systems based on their weak form is the fact that most of the common numerical methods in engineering, including commercial tools, are based on a Galerkin-type finite-dimensional approximation of the variational problem associated to the PDE. Assuming this perspective is even more natural for system classes, where the PDEs are derived from variational principles, as it is the case in structural mechanics 19 , see e.g. [START_REF] Reddy | Energy principles and variational methods in applied mechanics[END_REF].

Also in the context of existing works on linear port-Hamiltonian distributed-parameter systems in one spatial dimension, this perspective is natural. The statements on well-posedness and stability based on the theory of C 0 semigroups rely on the mild solution of the abstract (operator) differential equation. These solutions, however, corresponds to the weak solutions, as known from the theory of PDEs, see [START_REF] Jacob | Linear port-Hamiltonian systems on infinite-dimensional spaces[END_REF], page 127: "In fact, the concept of a mild solution is the same as the concept of a weak solution used in the study of partial differential equations."

A third point which motivates to discretize port-Hamiltonian distributed-parameter systems based on their weak form, is the close relation of discrete exterior calculus (i.e. the mathematical formalism for integral modeling of conservation laws), which has been used in [START_REF] Seslija | Explicit simplicial discretization of distributed-parameter port-hamiltonian systems[END_REF] for PH systems: "Note that the process of integration to suppress discontinuity is, in spirit, equivalent to the idea of weak form used in the Finite Element method." [START_REF] Mathieu Desbrun | Discrete differential forms for computational modeling[END_REF] 20 Finally, also in the work of Bossavit on the mixed geometric discretization for computational electromagnetism [START_REF] Bossavit | Differential forms and the computation of fields and forces in electromagnetism[END_REF], [START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements[END_REF], the quality of a weak formulation is addressed "How weak is the weak solution in finite element methods" [START_REF] Bossavit | How weak is the" weak solution" in finite element methods[END_REF].

In the following, we will lift the smoothness assumption of the differential forms and specify concrete functional spaces with conditions on their (weak) differentiability. This is necessary to choose appropriate basis forms for the finite-dimensional approximation of the Stokes-Dirac structure. For a quick overview on the necessary notions from functional analysis, we refer to Appendix ??.

The weak form of the Stokes-Dirac structure of Prop. 1 is obtained by a duality pairing (which involves exterior multiplication and integration) on Ω with test forms of appropriate degree which do not vanish on the boundary 21 . The latter allows for a weak imposition of the input boundary conditions e Γ i , i = 1, . . . , n Γ and êΓ j , j = 1, . . . , nΓ .

Definition 3. The weak form of the Stokes-Dirac structure of Prop. 1 is given by the subspace D ⊂ F ×E with 22

F = L 2 Λ p Ω × L 2 Λ q Ω × L 2 Λ n-p Γ1 × • • • × L 2 Λ n-p Γn Γ × L 2 Λ n-q Γ1 × • • • × L 2 Λ n-q Γn Γ E = H 1 Λ n-p Ω × H 1 Λ n-q Ω × L 2 Λ n-q Γ1 × • • • × L 2 Λ n-q Γn Γ × L 2 Λ n-p Γ1 × • • • × L 2 Λ n-p Γn Γ (31)
where

v p |f p Ω = v p |(-1) r de q Ω ∀v p ∈ H 1 Λ n-p (Ω), v q |f q Ω = v q |de p Ω ∀v q ∈ H 1 Λ n-q (Ω), (32)
and

f Γ i = tr Γi e p , e Γ i = (-1) p tr Γi e q , f Γ j = (-1) p tr Γj e q , êΓ j = tr Γj e p . ( 33
)
Remark 4. Note the subtle formal difference in the definition of the boundary port variables. In [START_REF] Fernandez-Nieto | Coupling superposed 1d and 2d shallow-water models: Source terms in finite volume schemes[END_REF] they are defined as restrictions of the smooth effort forms to the boundary, uniquely defined by continuity.

Here, we use the trace operator, which defines the extension of an H 1 form to the boundary, see Appendix ??. 19 By using a discrete variational formulation space-time discrete models can be even obtained directly, see e.g. [START_REF] Marsden | Multisymplectic geometry, variational integrators, and nonlinear pdes[END_REF]. 20 The paper is a good introduction to exterior differential calculus. 21 In the variational formulation of boundary value problems, mostly test functions with compact support inside Ω are chosen such that boundary conditions have to be imposed directly on the solution. This is however, not mandatory. By test functions which are non-zero on ∂Ω, boundary conditions can be imposed in a weak fashion, cf. [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF], Section 14.3.1, p. 483. 22 For brevity, the domain of the differential forms is written as an index, e.g.

L 2 Λ p Ω = L 2 Λ p Ω .
With the functional spaces of the flow and effort differential forms (f p , f q , f Γ 1 , . . . , f Γ n Γ , f Γ 1 , . . . , f Γ nΓ ) and (e p , e q , e Γ 1 , . . . , e Γ n Γ , êΓ 1 , . . . , êΓ nΓ ) integration by parts according to (??) is applicable, which yields the following weak form of the Stokes-Dirac structure.

Proposition 2. The weak form of the Stokes-Dirac structure in Prop. 1 with weak treatment of boundary port variables is given by the subset D ⊂ F × E, F and E according to [START_REF] Hecht | New development in freefem++[END_REF], where

v p |f p Ω = (-1) r+q dv p |e q Ω -(-1) r+p+q n Γ i=1 tr v p |e Γ i Γi -(-1) r+p+q nΓ j=1 tr v p | f Γ j Γj v q |f q Ω = (-1) p dv q |e p Ω - (-1) p n Γ i=1 tr v q |f Γ i Γi - (-1) p nΓ j=1 tr v q |ê Γ j Γj , ( 34
)
holds for all test forms v p ∈ H 1 Λ n-p (Ω) and v q ∈ H 1 Λ n-q (Ω).

Proof. Eq. ( 34) follows from [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF] via integration by parts and the identities

(-1) p v p |tr e q ∂Ω = n Γ i=1 v p |e Γ i Γi + nΓ j=1 v p | f Γ j Γj , v q |tr e p ∂Ω = n Γ i=1 v q |f Γ i Γi + nΓ j=1 v q |ê Γ j Γj . ( 35
)
The latter are due to the definition (33) of boundary port variables and the definition of the subsets Γ i , Γj which cover completely ∂Ω and whose intersections are either empty or have dimension n -1.

Remark 5. The latter representation of the Stokes-Dirac structure -if considered on a single control volume -is suitable for discontinuous Galerkin schemes, see e.g. [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF], where the boundary terms are replaced by suitable numerical fluxes.

Using the effort forms as test forms, v p = e p , v q = e q , and adding both equations of (34), we obtain (after some reformulations and exploiting [START_REF] Geuzaine | GetDP: a general finite-element solver for the de Rham complex[END_REF], see appendix) e p |f p Ω + e q |f q Ω + (-1) p e q |e p ∂Ω = 0, [START_REF] Kato | Perturbation theory for linear operators[END_REF] or, with the definition of boundary port variables,

e p |f p Ω + e q |f q Ω + n Γ µ=1 e ∂,µ |f ∂,µ Γµ + nΓ ν=1 f ∂,ν |ê ∂,ν Γν = 0, (37)
which corresponds to the initially derived power continuity equation [START_REF] Luiz Cardoso-Ribeiro | A port-Hamiltonian model of liquid sloshing in moving containers and application to a fluid-structure system[END_REF].

We have arrived at a weak (or variational) representation of the Stokes-Dirac structure of Prop. (1) which suits to establish discretized mixed Galerkin models of PH systems of two conservation laws.

Mixed Galerkin approximation

In this section, we introduce the mixed Galerkin approximation of the weak form of the Stokes-Dirac structure for a system of two conservation laws. For convenience of notation, we omit to explicitly write out the trace operator on the subsets of the boundary, i.e. v p |e q Γi := tr v p |tr e q Γi etc. in the sequel. We start with the representation

v p |f p Ω = (-1) r+q dv p |e q Ω -(-1) r+q n Γ i=1 v p |e q Γi -(-1) r+q nΓ j=1 v p |e q Γj v q |f q Ω = (-1) p dv q |e p Ω -(-1) p n Γ i=1 v q |e p Γi -(-1) p nΓ j=1 v q |e p Γj , (38)
i.e. [START_REF] Holm | Geometric mechanics: Dynamics and symmetry[END_REF] without the explicit denomination of the boundary port variables. For a mixed Galerkin approximation of the Stokes-Dirac structure, we

• use different (mixed) bases to approximate the spaces of flow and effort forms

• and, from these bases, we choose the appropriate ones to approximate the test forms (Galerkin method).

For the Stokes-Dirac structure, the effort bases are the natural choice to approximate the spaces of test forms, as they retain the interpretation of the duality products as power products, see [START_REF] Kato | Perturbation theory for linear operators[END_REF].

Note that from the original (strong) formulation of the Stokes-Dirac structure, until its weak formulation, both equations feature only the exterior derivative and not the co-differential, i.e. no information about the metric on the functional spaces is contained. Nor do we decide -for a finite-element (FE) discretization which serves as illustrative example -to approximate the two equations on different, staggered grids, as it is common practice in computational electromagnetism [START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements[END_REF] or the numerical treatment of transport phenomena [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]. Instead, we approximate the variables in both conservation laws with the same conforming subspaces (i.e. on the same mesh in FE), which has the advantage that boundary variables are defined directly on ∂Ω, without having to cope with an eventual grid shift.

Nevertheless, the power-preserving maps that will be introduced to ensure non-degenerate duality pairings between the discrete bond variables will give rise to an interpretation of the resulting schemes in terms of staggered grids by the localization of the mapped discrete variables. The ability to parametrize the power-preserving maps will allow to adapt the approximation to the nature of the system, e.g. to account for diffusive effects in a hyperbolic system.

Approximation problem and compatibility condition

The flow differential forms will be approximated by linear combinations of appropriate basis forms of subspaces

Ψ p h = span{ψ p 1 , . . . , ψ p Np } ⊂ L 2 Λ p (Ω) Ψ q h = span{ψ p 1 , . . . , ψ q Nq } ⊂ L 2 Λ q (Ω). (39)
The subspaces for the effort and test forms are, accordingly,

Φ p = span{ϕ p 1 , . . . , ϕ p Mp } ⊂ H 1 Λ n-p (Ω) Φ q h = span{ϕ q 1 , . . . , ϕ q Mq } ⊂ H 1 Λ n-q (Ω). ( 40
)
From the trace theorem for H 1 spaces (see Appendix ??), we know that the restriction of the latter spaces to the boundary is L 2 . h > 0 denotes the discretization parameter 23 and we assume an appropriate choice of approximation spaces, i.e. for a given functional space V and its approximation V h (see [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF], Section 5.2) it is true that inf

v h ∈V h v -v h → 0 for all v ∈ V if h → 0.
The mixed Galerkin approximation problem is as follows: Find approximate flow and effort forms

f p h (z) = Np k=1 f p k ψ p k (z) = f p , ψ p (z) ∈ Ψ p h , f q h (z) = Nq l=1 f q l ψ q l (z) = f q , ψ q (z) ∈ Ψ q h , (41)
and

e p h (z) = Mp i=1 e p i ϕ p i (z) = e p , ϕ p (z) ∈ Φ p h , e q h (z) = Mq j=1
e q j ϕ q j (z) = e q , ϕ q (z) ∈ Φ q h , (42) 23 Which corresponds to the spatial extent of finite elements or the reciprocal order of a polynomial approximation.

where •, • denotes the standard inner product in R n , such that

v p h |f p h Ω = (-1) r+q dv p h |e q h Ω -(-1) r+q n Γ µ=1 v p h |e q h Γµ -(-1) r+q nΓ ν=1 v p h |e q h Γν v q h |f q h Ω = (-1) p dv q h |e p h Ω -(-1) p n Γ µ=1 v q h |e p h Γµ -(-1) p nΓ ν=1 v q h |e p h Γν . ( 43
)
holds for all v p h ∈ Φ p h , v q h ∈ Φ q h . The discrete flow and effort vectors

f p = [f p 1 , . . . , f p Np ] T , f q = [f q 1 , . . . , f q Nq ] T ,
and e p = [e p 1 , . . . , e p Mp ] T , e q = [e q 1 , . . . , e q Mq ] T , [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] contain the approximation coefficients, and the vectors (we omit the argument z in the sequel)

ψ p (z) = [ψ p 1 (z), . . . , ψ p Np (z)] T , ψ q (z) = [ψ q 1 (z), . . . , ψ q Nq (z)] T ,
and

ϕ p (z) = [ϕ p 1 (z), . . . , ϕ p Mp (z)] T , ϕ q (z) = [ϕ q 1 (z), . . . , ϕ q Mq (z)] T , ( 45
)
contain the approximation basis forms.

The flow variables are understood as time derivatives of the distributed states with negative sign, see [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]. Thus, they are approximated with the same bases,

p h (z) = Np k=1 p k ψ p k (z) = p, ψ p (z) ∈ Ψ p h , q h (z) = Nq l=1 q l ψ q l (z) = q, ψ q (z) ∈ Ψ q h , (46)
and p = [p 1 , . . . , p Np ] T and q = [q 1 , . . . , q Nq ] T . (

denote the discrete state vectors.

The mixed Galerkin approximation (43) of ( 38) is exact for flow and effort forms from the approximation spaces (39), [START_REF] Macchelli | Energy shaping of distributed parameter port-hamiltonian systems based on finite element approximation[END_REF] (in these subspaces, the residual error vanishes), if the following compatibility conditions hold:

span{ψ p 1 , . . . , ψ p Np } = span{dϕ q 1 , . . . , dϕ q Mq } span{ψ q 1 , . . . , ψ q Nq } = span{dϕ p 1 , . . . , dϕ p Mp }.

In contrast to [START_REF] Golo | Hamiltonian discretization of boundary control systems[END_REF] (Assumptions 3 and 7), this compatibility of forms is understood in the weak sense. This means, more precisely -consider the original weak formulation [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF] and the definition of the weak exterior derivative -that for all test forms with compact support inside Ω, i.e.

v p ∈ H 1 0 Λ n-p (Ω), v q ∈ H 1 0 Λ n-q (Ω), there exist constants a p k , a q l , b p i , b q j such that Np k=1 a p k v p |ψ p k Ω + Mq j=1 b q j v p |dϕ q j Ω = 0, Nq l=1 a q l v q |ψ q l Ω + Mp i=1 b p i v q |dϕ p i Ω = 0. ( 49
)

Approximation of the variational problem

The weak formulation [START_REF] Kotyczka | Discrete port-Hamiltonian formulation and numerical approximation for systems of two conservation laws[END_REF] of the Stokes-Dirac structure represents a variational problem: flow and effort forms are sought under arbitrary variation of the test forms. Replacing the finite-dimensional approximations ( 41), [START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF] and using the test forms from the effort bases,

v p h = v p , ϕ p , v q h = v q , ϕ q , v p ∈ R Mp , v q ∈ R Mq , ( 50
)
we obtain the equations (the exterior derivative applies element-wise to a vector of differential forms)

v p , ϕ p f p , ψ p Ω -(-1) r+q v p , dϕ p e q , ϕ q Ω + (-1) r+q n Γ µ=1 v p , ϕ p e q , ϕ q Γµ + (-1) r+q nΓ ν=1
v p , ϕ p e q , ϕ q Γν = 0 (51) and v q , ϕ q f q , ψ q Ω -(-1) p v q , dϕ q e p , ϕ p

Ω + (-1) p n Γ µ=1 v q , ϕ q e p , ϕ p Γµ + (-1) p nΓ ν=1 v q , ϕ q e p , ϕ p Γν = 0 (52)
Evaluating the integrals over the basis forms, we get the short notation of the finite dimensional variatonal problem

v p , M p f p + v p , K p + n Γ µ=1 L µ p + nΓ ν=1 Lν p e q = 0 ∀v p ∈ R Mp , v q , M q f q + v q , K q + n Γ µ=1 L µ q + nΓ ν=1
Lν q e q = 0 ∀v q ∈ R Mq .

(

) 53
whose matrices

M p ∈ R Mp×Np , M q ∈ R Mq×Nq , K p , L µ p , Lν p ∈ R Mp×Mq , K q , L µ q , Lν q ∈ R Mq×Mp , µ = 1, . . . , n Γ , ν = 1, . . . , nΓ , are composed of the elements [M p ] ik = ϕ p i |ψ p k Ω , [M q ] jl = ϕ q j |ψ q l Ω , [K p ] ij = -(-1) r+q dϕ p i |ϕ q j Ω , [K q ] ji = -(-1) p dϕ q j |ϕ p i Ω , [L µ p ] ij = (-1) r+q ϕ p i |ϕ q j Γµ , [L µ q ] ji = (-1) p ϕ q j |ϕ p i Γµ , [ Lν p ] ij = (-1) r+q ϕ p i |ϕ q j Γν , [ Lν q ] ji = (-1) p ϕ q j |ϕ p i Γν . ( 54
)
Equations ( 53) must hold for arbitrary vectors (discrete variations) v p , v q , which gives the matrix equations

M p f p + (K p + L p )e q = 0, M q f q + (K q + L q )e p = 0. (55)
It is straightforward to show by application of (??) that

[L µ p ] ij = [L µ q ] ji , [ Lν p ] ij = [ Lµ q ] ji (56)
i.e. L µ p = (L µ q ) T and Lν p = ( Lν q ) T . By defining

L p = n Γ µ=1 L µ p + nΓ ν=1 Lν p , L q = n Γ µ=1 L µ q + nΓ ν=1 Lν q , ( 57
)
we can prove the following.

Proposition 3.

The following relation holds between the elements of the matrices K p , K q and L p , L q .

[K p + L p ] ij + [K q + L q ] ji = [L p ] ij = [L q ] ji , (58)
i.e. (K p + L p ) + (K q + L q ) T = L p = L T q .

Proof. By the definition (57), the elements of L p , L q are duality products over the effort basis forms on the complete boundary ∂Ω. Thus, we have that

[K p + L p ] ij + [K q + L q ] ji = -(-1
) r+q dϕ p i |ϕ q j Ω + (-1) r+q ϕ p i |ϕ q j ∂Ω -(-1) p dϕ q j |ϕ p i Ω + (-1) p ϕ q j |ϕ p i ∂Ω . (59) Using skew-symmetry of the wedge product and the integration-by-parts formula for differential forms, the right hand side can be rewritten as

(-1) p ϕ q j (z)|ϕ p i (z) ∂Ω = [L q ] ji = [L p ] ij , ( 60
)
which completes the proof.

Proposition 4. The finite-dimensional approximation of the power continuity equation ( 36) is e p , M p f p + e q , M q f q + e q , L q e p = 0 (61a) or e p , M p f p + e q , M q f q + L p e q , e p = 0, (61b

)
where •, • denotes the standard inner product on Euclidean space.

Proof. The relations follow from replacing v p = e p , v q = e q in the finite-dimensional variational problem [START_REF] Schoberl | Analysis and comparison of port-hamiltonian formulations for field theories-demonstrated by means of the mindlin plate[END_REF], adding both equations, substituting the definition (57) of L p , L q , and exploiting the property [START_REF] Trenchant | Structure preserving spatial discretization of 2D hyperbolic systems using staggered grids finite difference[END_REF].

As the matrices L p , L µ p , Lν p and L q , L µ q , Lν q describe the approximate power transferred over the boundary ∂Ω or parts of it, we refer to them as them boundary power matrices. L p = L T q , whose elements are the integrals of duality pairings between approximation forms on the boundary, will have reduced rank. The reason is that basis forms for interior effort degrees of freedom will be, in general, zero on the boundary24 .

Discrete boundary port variables

In the next step, we characterize mappings to discrete boundary port variables, according to the subdivision of the boundary ∂Ω into disjoint regions Γ i , i = 1, . . . , n Γ , Γj , j = 1, . . . , nΓ . These subsets of the discrete effort variables will be assigned either the role of discrete boundary in-or outputs. , associated with the boundary subdomains Γ µ ⊂ ∂Ω, µ = 1, . . . , n Γ , Γν ⊂ ∂Ω, ν = 1, . . . , nΓ are defined via the discrete expression of the power transfered over Γ µ and Γν , respectively: e q , L µ q e p =: e b,µ , f b,µ 0 , e p , Lν p e q =: êb,ν , f b,ν

0 . ( 62
)
We decompose the boundary power matrices for each boundary subdomain in matrix products

L µ q = (T µ q ) T S µ p,0 , Lν p = ( Tν p ) T Ŝν q,0 . ( 63
)
The boundary trace matrices

26 T µ q ∈ R M µ b ×Mq , Tν p ∈ R M ν b ×Mp define the effort degrees of freedom e b,µ = T µ q e q , êb,ν = Tν p e p . ( 64
)
that lie on the boundary and are assigned the roles of input variables. We call S µ p,0 ∈ R M µ b ×Mp , Ŝν q,0 ∈ R M ν b ×Mq the collocated boundary output matrices. They define the boundary flow variables

f b,µ 0 = S µ p,0 e p , f b,ν 0 = Ŝν q,0 e q , ( 65
)
which, together with the discrete efforts [START_REF] Minh | Symplectic spatial integration schemes for systems of balance equations[END_REF], satisfy exactly the discrete power balance (62) on the different portions of the boundary.

Because of e ∂,µ |f ∂,µ Γµ = (-1) p e q |e p Γµ ≈ (-1) p e q , ϕ q e p , ϕ p Γµ = e q , L µ q e p = e b,µ , f b,µ 0 , ê∂,ν | f ∂,ν Γν = (-1) p e p |e q Γν ≈ (-1) p e p , ϕ p e q , ϕ q Γν = e p , Lν p e q = êb,ν , f b,ν 0 , [START_REF] Whitney | Geometric integration theory[END_REF] the definition of discrete boundary port variables is consistent with the distributed definition [START_REF] Rr Hiemstra | High order geometric methods with exact conservation properties[END_REF]. Summation over the individual boundary power matrices according to [START_REF] Tonti | A direct discrete formulation of field laws: The cell method[END_REF], yields the boundary power balance in compact notation

T T q S p,0 + ŜT q,0 Tp = L q , ( 67
)
where

T q =    T 1 q . . . T n Γ q    , S p,0 =    S 1 p,0 . . . S n Γ p,0    , Ŝq,0 =    Ŝ1 q,0 . . . Ŝn Γ q,0    , Tp =    T1 p . . . Tn Γ p    . ( 68
)
The discrete power continuity equation, which represents the counterpart of ( 37) in the approximation subspaces, finally reads e p , M p f p + e q , M q f q + e b ,

f b 0 + êb , f b 0 = 0, ( 69
)
where the vectors

e b =    e b,1 . . . e b,n Γ    , f b 0 =    f b,1 0 . . . f b,n Γ 0    , êb =    e b,1 . . . e b,n Γ    , f b 0 =    f b,1 0 . . . f b,n Γ 0    . ( 70
)
satisfy

e b = T q e q , f b 0 = S p,0 e p , êb = Tp e p , f b 0 = Ŝq,0 e q . ( 71
)
We are at the point to express the subspace of the discrete bond space on which power continuity according to (69) holds, corresponding to the infinite-dimensional Stokes-Dirac structure in Prop. 1:

Proposition 5. Given the finite-dimensional space of discrete flows and efforts 53) and ( 64), ( 65) hold} (73)

F × E = R Np × R Nq × R M b × R Mb × R Mp × R Mq × R M b × R Mb (72) with M b = n Γ µ=1 M µ b , Mb = nΓ ν=1 M ν b . The subspace {(f p , f q , f b 0 , f b 0 , e p , e q , e b , êb ) ∈ F × E | (
satisfies the isotropy 27 condition, D ⊂ D ⊥ with respect to the bilinear form •, • M that results from symmetrization of the degenerate duality product

•, • M := e p , M p f p + e q , M q f q + e b , f b 0 + êb , f b 0 . ( 74
)
Proof. Along the lines of Prop. 18 in [START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF].

The subspace described in the proposition is, however, not a Dirac structure. This is due to the fact that the duality product •, • M is degenerate in general, i.e. its value can be zero for nonzero discrete flows and/or efforts that lie in the kernel of M p , M q , or their transposes 28 .

27 A quadratic form is called isotropic if there exists a non-zero vector which is mapped to zero. 28 Consider systems of conservation laws with p = n, q = 1. In 1D, we have for Whitney finite elements on a simplicial triangulation or a polynomial approximation [START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF], dim e p = dim f p + 1 and dim e p = dim f p + 1, i.e. dim(ker M T p ) = 1 and dim(ker M T q ) = 1. On a 2D triangulation with Whitney forms, Mq will be square, while Mp will be fat, for a certain minimal grid size, i.e. ker Mp = ∅.

Discrete conservation laws

Instead of the mixed Galerkin approximation of (32) on the complete spatial domain Ω, we can approximate both conservation laws on the discrete, oriented geometric objects of dimension p and q (faces and edges in 2D, for example) of the discretization mesh. We obtain a discrete representation of the two conservation laws in the form

f p f q = 0 (-1) r d p d q 0 e p e q . ( 75
)
D p and D q are considered discrete derivative matrices, which replace the exterior derivative in the original infinite-dimensional representation. If Whitney forms of lowest polynomial degree are used as finite element approximation bases, they represent the co-incidence or co-boundary relations between the discrete geometric objects associated to effort and flow forms. The term co-incidence matrix is preferred to co-boundary matrix in the PH framework to avoid confusion with the system boundary, see [START_REF] Seslija | Explicit simplicial discretization of distributed-parameter port-hamiltonian systems[END_REF] or [START_REF] Van Der Schaft | Discrete conservation laws and port-hamiltonian systems on graphs and complexes[END_REF].

Replacing (75) in [START_REF] Seslija | Explicit simplicial discretization of distributed-parameter port-hamiltonian systems[END_REF], we obtain the relations

(-1) r M p d p + (K p + L p ) = 0 M q d q + (K q + L q ) = 0 (76)
for the matrices of both representations.

Geometric discretization of the port-Hamiltonian system

By power-preserving projections of the discrete power variables (flows and efforts), we formulate an approximate Dirac structure, which is based on non-degenerate duality pairings between discrete bond variables. This Dirac structure admits different representations, in particular it can be written in an explicit input-output form. The consistent approximation of the constitutive relations, which defines a discrete Hamiltonian, finally yields the family of PH finite-dimensional state space models, which are parametrized by the power-preserving projections and the definitions of boundary control inputs.

Power-preserving projections and conjugated output maps

Due to the different approximation bases for flows and efforts, the matrices M p and M q can be nonsquare such that the duality products e p , M p f p and e q , M q f q are degenerate. An approximation Dirac structure must be based on a non-degenerate duality product between the variables of the discrete bond space. To this end, we will determine power-preserving projections ẽp = P ep e p , ẽq = P eq e q and f p = P f p f p , f q = P f q f q .

(77) such that

Ñp := dim ẽp = dim f p ≤ rank(M p ) and Ñq := dim ẽq = dim f q ≤ rank(M q ). ( 78
)
We refer to the vectors f p , ẽp ∈ R Ñp , f q , ẽq ∈ R Ñq as projected discrete flows and efforts, as they can be interpreted as projections on subspaces of the original discrete bond spaces.

Example 3. In the 1D case, p = q = 1, using Whitney finite elements or the pseudo-spectral method [START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF], we have, N = N p = N q and M = M p = M q with M = N + 1. Fixing f p = f p , f q = f q and P ep = M T p , P eq = M T q , the new efforts ẽp , ẽq are projections of e p , e q in direction of the kernel of M T p and M T q , respectively.

The following definition summarizes the core property of power-preserving projections.

Definition 5. The discrete flow and effort projections (77) are power-preserving if they satisfy a discrete power balance ẽp ,

f p + ẽq , f q + e b , f b + êb , f b = 0 (79)
with the given boundary inputs e b , êb according to [START_REF] Minh | Symplectic spatial integration schemes for systems of balance equations[END_REF] and possibly modified boundary outputs

f b = S µ p e p , f b = Ŝν q e q . ( 80
)
Remark 6. If the projections satisfy P T ep P f p = M p and P T eq P f q = M q , the "interior" part of the power balance ( 61) is exactly represented by the new flows and efforts, and (79) holds with the original, collocated outputs f b 0 , f b 0 . If, however, Ñq < rank(M q ) and/or Ñp < rank(M p ), a part of the power, originally described by e p , M p f p + e q , M q f q , must be "swapped" to the boundary terms of (79) via the redefinition of the outputs. This way, the power-balance is maintained globally, and conservativeness of the finite-dimensional approximation is guaranteed.

To characterize the power-preserving projections and modified output maps that guarantee power continuity (79), we replace the definitions of projections, in-and outputs, and substitute f p , f q according to the discrete representation (75) of the conservation laws. The new power variables are now expressed in terms of the original discrete efforts,

    f p f q f b f b     f =     0 (-1) r P f p d p P f q d q 0 0 Ŝq S p 0     E T e p e q e ,     ẽp ẽq êb e b     ē =     P ep 0 0 P eq Tp 0 0 T q     F T e p e q e . (81)
Equation ( 79) must hold for arbitrary e p , e q , and we obtain the following matrix condition.

Proposition 6. The effort and flow projections and output maps are power-preserving, if they satisfy the matrix equation (-1) r d T p P T f p P ep + P T eq P f q d q + T T q S p + ŜT q Tp = 0.

The power-preserving projections and output maps are not unique. Different parametrization of the matrices, yield different finite-dimensional Dirac structures that approximate the original Stokes-Dirac structure of Proposition 1. Together with a consistent approximation of the constitutive equations, we can obtain PH approximate models with different numerical properties. A favorable parametrization will depend on the nature of the system (e.g. if the closure equations make the system hyperbolic or parabolic), the distribution and type of boundary inputs, and the application case.

In Section 5, we will illustrate the power-preserving projections at the example of Whitney approximation forms on a rectangular simplicial mesh in 2D. The degrees of freedom in the projections will allow for a trade-off between centered schemes and upwinding in the discretized PH models.

Dirac structure on the new bond space

The power-preserving projections and output maps that satisfy (82) define a Dirac structure. We verify that (81) is an image representations of this Dirac structure on the new discrete bond space. If the effort maps are invertible, an unconstrained input-output representation exists.

Proposition 7 (Image representation). Consider the discrete flow and effort vectors f and ē as indicated in (81). ( f , ē) is an element of the bond space F × Ē = R Ñp+ Ñq+Mb+ Mb × R Ñp+ Ñq+Mb+ Mb . Let Ñp + Ñq + Mb + M b = M p + M q and the projection and output matrices satisfy the matrix condition (82). If

rank( P ep Tp ) = M p and rank( P eq T q ) = M q , ( 83
)
then the subspace D = {( f , ē) ∈ F × Ē | f = E T e, ē = F T e, e ∈ R M p +M q }, (84)
is a Dirac structure.

Proof. According to the definition of the image representation of a Dirac structure (see e.g. [START_REF] Van Der Schaft | L2-Gain and Passivity Techniques in Nonlinear Control[END_REF], Section 4.4.1) the dimensions of f and ē must be less 29 or equal dim(e), which is ensured by Ñp + Ñq + Mb + M b = M p + M q . The condition rank([F E]) = M p + M q is satisfied by (83), from which rank(F) = M p + M q follows. Moreover, the skew-symmetry condition EF T + FE T = 0 must hold. EF T + FE T according to (81) gives 0 d T q P T f q P eq + (-1) r P T ep P f p d p + S T p T q + ŜT p Tq (-1) r d T p P T f p P ep + P T eq P f q d q + T T q S p + ŜT q Tp 0,

which equals zero by the characterization of power-preserving projections and output maps (82).

Corollary 2 (Input-output representation). Under the conditions of Prop. 7, the Dirac structure admits an unconstrainted input-output representation

-f p -f q = 0 J p J q 0 J=-J T ẽp ẽq + 0 B p B q 0 B êb e b f b f b = 0 C q C p 0 C=B T ẽp ẽq + 0 D q D p 0 D=-D T êb e b . ( 86
)
with

J p = -J T q , B p = C T p , B q = C T q , D q = D T p . ( 87
)
Proof. The (skew-)symmetry conditions can be summarized as

-J p -B p C q D q + -J q -B q C p D p T = 0 0 0 0 . ( 88
)
By invertibility of the matrices in (83),

-J p -B p C q D q = (-1) r P f p d p Ŝq P eq T q -1 , -J q -B q C p D p = P f q d q S p P ep Tp -1 . ( 89
)
Substituting these relations in (88) and multiplying with P T ep TT p from the left and P eq T q from the right, yields the left hand side of (82). The right hand side being zero, this proves (skew-)symmetry of the matrices (87) of the input-output representation.

The proposition is a generalization of Prop. 20 in [START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF] for the 1D case and the pseudo-spectral method. Note that the rank condition (83) on the effort and flow and boundary maps is sufficient (not necessary) for the subspace (84) to be a Dirac structure. The fact that both matrices in (83) are assumed square and invertible, guarantees the input-output representation in the corollary.

Finite-dimensional port-Hamiltonian model

To render the input-output representation of the Dirac structure a finite-dimensional port-Hamiltonian model of a canonical system of two conservation laws, we replace the projected discrete flow variables by time derivatives of discrete states 30-

f p =: ṗ ∈ R Ñp , -f q =: q ∈ R Ñq . ( 90
)
On the other hand, the projected efforts need to be replaced by the partial derivatives of a suitable discrete Hamiltonian Hd (p, q)

ẽp = ∂ Hd ∂ p T ∈ R Ñp , ẽq = ∂ Hd ∂ q T ∈ R Ñq . ( 91
)
The discrete Hamiltonian must be defined in such a way that the discrete effort variables represent a consistent approximation of their continuous counterpart. We present the discretization of the constitutive equations in detail in the FE example of Section 5.

With the state, input and output vectors

x = p q , u = êb e b , y = f b f b , ( 92
which is the finite-dimensional counterpart of [START_REF] Baaiu | Structure-preserving infinite dimensional model reduction: Application to adsorption processes[END_REF]. The PH form allows to easily interconnect the finitedimensional model of the system of two conservation laws with other subsystems in a power-preserving way, which is the basis for energy-based control design by interconnection see e.g. [START_REF] Macchelli | Energy shaping of distributed parameter port-hamiltonian systems based on finite element approximation[END_REF].

Can we use the projections to define new approximation bases that can be used in the variational problem and directly yield square and invertible matrices M p , M q ?

Example: Whitney forms on a 2D rectangular grid

To illustrate the steps to obtain an approximate PH state space model with desired boundary inputs by geometric discretization, we consider a 2-dimensional rectangular domain Ω = (0, L x ) × (0, L y ) ⊂ R 2 , with boundary ∂Ω, covered by a regular, oriented simplicial triangulation T h , as sketched in Fig. 2. The approximation bases for flows and efforts [START_REF] Le Gorrec | Dirac structures and boundary control systems associated with skew-symmetric differential operators[END_REF], [START_REF] Macchelli | Energy shaping of distributed parameter port-hamiltonian systems based on finite element approximation[END_REF] are composed of Whitney forms [START_REF] Whitney | Geometric integration theory[END_REF] of lowest polynomial degree, which can be constructed based on the barycentric node weights [START_REF] Bossavit | Generating whitney forms of polynomial degree one and higher[END_REF]. The degrees of freedom are directly associated to the nodes, directed edges and faces of the mesh. The well-known geometric discretization of Maxwell's equations [START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements[END_REF] is based on Whitney forms, and the resulting finitedimensional models feature the (co-)incidence matrices of the underlying discretization meshes [START_REF] Bossavit | Yee-like schemes on a tetrahedral mesh, with diagonal lumping[END_REF]. They can be considered a direct representation of the physical laws on the discrete balance regions of the triangulation, which underscores the strong relation to Tonti's cell method [START_REF] Tonti | A direct discrete formulation of field laws: The cell method[END_REF]. In contrast to [START_REF] Bossavit | Yee-like schemes on a tetrahedral mesh, with diagonal lumping[END_REF], [START_REF] Tonti | A direct discrete formulation of field laws: The cell method[END_REF], where the conservation laws are evaluated on dual or staggered grids, we start with a single mesh. Nevertheless, in our approach, the projections of the original degrees of freedom gives rise to interpret the projected flows and efforts in terms of topological duality.

Mesh, matrices and dimensions

Using Whitney basis forms, the degrees of freedom in the mixed Galerkin approach are associated to integrals of distributed quantities on the k-simplices of the mesh. The dimensions of the (initial) discrete flow and effort vectors equal the numbers of corresponding nodes, edges and faces on the grid. The same holds for the discrete efforts on the boundary, which are designated inputs and are localized at the corresponding boundary nodes and edges, see Table 1.

The mixed Galerkin approximation of the Stokes-Dirac structure yields a set of matrices with different sizes and ranks, see Table 2. The construction of power-preserving projections and conjugated output matrices that satisfy the matrix equation ( 82), is based on rank considerations of the involved matrix products.

Power-preserving projections, in-and outputs

We describe the different matrices that define the discrete power variables on the projected bond space.

Input trace matrices and effort projections

Identifying the elements of the input vector u = (ê b ) T (e b ) T T with effort degrees of freedom on the boundary nodes and edges corresponds to a consistent imposition of the effort boundary conditions in the finite-dimensional model. To arrive at the input-output representation (86), the matrices Π p := P ep Tp and Π q := P eq T q (95) should be square and invertible. Choosing Π p and Π q as permutation matrices, each projected effort/boundary effort is assigned to an interior/boundary simplex. The property Π -1 p/q = Π T p/q makes the matrices of the state space model, which result from a right multiplication with Π -1 p/q , particularly simple.

Flow/state projections

For the considered rectangular simplicial grid, we propose the following constructions for P f p and P f q . We omit the derivations, and illustrate the results graphically.

Each element f p i , i = 1, . . . , Ñp , of f p = P f p p is related to a 2-chain (a weighted formal sum of 2simplices), located around the node associated to ẽp i . The node and the weighted 2-chain can be considered, hence, topologically dual objects. The weights of the 2-simplices around a node, are illustrated in Fig. 3. The index I refers to the lower right, II to the upper left triangles and

α ν + β ν + γ ν = 1, (j -1)(N +1) +k (j -1)(N +1) +k+1 j(N +1) +k j(N +1) +k+1 (j -1)N +k jN +k (j -1)(N +1) +K +k (j -1)(N +1) +K +k+1 (j -1)N +L + k (j -1)N +k (M +j -1)N +k 1 2 3
) 96
where p is a weighted sum of the original discrete states p j on 2-simplices around boundary input nodes.

p describes a defect of the total conserved quantity on the original grid. This defect at the boundary must be considered in the consistent definition of the discrete constitutive equations as shown further below.

Remark 7. By the topological duality between the elements of f p and ẽp as illustrated in Fig. 3, it is easy to imagine that the deviation p can be related to well-known effects from the discretization with staggered grids, like ghost values, see e.g. [START_REF] Kotyczka | Discrete port-Hamiltonian formulation and numerical approximation for systems of two conservation laws[END_REF] for a discussion from the PH point of view.

A related interpretation holds for the different elements of the projected flow vector f q = P f q f q . As shown in Fig. 4, each element of f q can be considered dual to a discrete effort ẽq i on a horizontal, vertical or diagonal edge (drawn in red). f p i is constructed as a linear combination of the elements of f q on adjacent edges. The weights can be clustered, which corresponds to a decomposition of f q i into a rotational part and components along and across the effort edge. Only the latter part (drawn in blue) contributes to the discrete constitutive equations as discussed further below.

Conjugated discrete outputs

Like the projected flows and efforts, the discrete outputs f b = S p e p and f b = Ŝq e q are constructed as weighted sums of the discrete efforts in the vicinity of the dual boundary input, see Fig. 5. f b i is defined by a convex sum of node efforts. Note that f b i is composed of a "rotation" part and a convex sum associated to the direct neighbor edges. Remark 8. By reconstructing the rotational components of q from the given quantities, can be used to discretize the vorticity term in the shallow water equation [START_REF] Bossavit | Generating whitney forms of polynomial degree one and higher[END_REF].

Show by short calculation that the models indeed have no feedthrough.

Discrete constitutive equations

To obtain a consistent numerical approximation of the system of conservation laws, the discrete states p, q and the above-defined efforts ẽp , ẽq must be related via discrete constitutive relations that are consistent with the continuous ones. We consider the case of linear constitutive equations e p = δ p H = * p, e q = δ q H = - * q.

(97) that are derived from the simple quadratic Hamiltonian

ε II 3(εI +εII ) -(βI+βII ) ε I ε II -β II -ε II -ε I ε I -β I = ε I ε I ε I ε I ε II ε II ε II ε II + ẽq,a i 1-αI -αII -β II -β I (a) Weights for f q,a i -γI 2 -γII 2 -γII 2 -γI 2 ẽq,c k ≈ - γ I - γ I I (b) Weights for f q,c i -δ I +α I δ II δ I -δ II +α II -δ I -3(δI +δII ) +αI +αII -δ II = δ I δ I δ I δ I δ II δ II δ II δ II + -α I -α II
H = Ω ( 1 2 p ∧ * p + 1 2 q ∧ * q). (98)
The discrete constitutive equations will be expressed by ẽp = Q p p, ẽq = Q q q (99) with positive definite, diagonal matrices Q p , Q q . p and q are constructed as linear combinations of integral conserved quantities on the 2-and 1-simplices of the discretization grid. By the construction illustrated above, we can think p, q localized on a (virtual) dual grid, whose localization is parameterized by the projection parameters α I/II and β I/II . This justifies the interpretation of Q p , Q q as (diagonal)

discrete Hodge matrices that replace the Hodge star in (97).

For the consistent construction of these Hodge matrices, we assume a steady state. In this case, the elements of ẽp must represent "average" values of p on the weighted balance areas 32 for the pi . Accordingly, the elements of ẽq must reflect the integral flux of the vector field 33 q across the corresponding horizontal, vertical or diagonal edges. The weights in the definition of qi that belong to edges perpendicular to the ẽq i determine the diagonal elements of Q -1 q .

Simulation study

We consider the linear wave equation with Hamiltonian density (98) on a square domain Ω = [0, 20]×[0, 20] to illustrate the effects of different projection parameters, see Table 3. We impose the boundary conditions e p (0, 0) = u(t) = sin 2 ( π 8 t), 0 ≤ t < 8 0 t ≥ 8 , e q (x, y) = 0 on ∂Ω. Fig. 6 shows the simulated propagation of the wave in radial direction under different parametrizations of the method, see Table 3. The red line displays a circle with radius t sim -T /2 = 16, as a reference for the maximum of the wavefront 34 . The first parameter set produces a centered numerical scheme with respect to the propagation direction. The three remaining parameter sets correspond to a stronger weighting of the simplices in the propagation direction for the computation of f p i , see 3. As a result, the numerical schemes can be considered different realizations of "upwinding". Observe also the different variation from the shape of a quarter circle, which are due to the non-isotropic mesh and the different weights α I/II , β I/II .

Spectral approximation

For a short analysis of the spectral approximation properties, we consider the 1D wave equation on a domain Ω = [0, 1], which is characterized by the degree p = 1 (instead of p = 2) of the first conserved quantity differential form. With the quadratic Hamiltonian as in (98), the sign of the second constitutive equation e q = * q changes. The exponent in the canonical differential operator now becomes r = pq+1 = 2. Table 4 shows the approximation of the spectrum of the canonical differential operator of the Stokes-Dirac structure (exact eigenvalues ± 2k-1 2 πi, k = 1, 2, 3, . . .) under homogeneous Dirichlet boundary conditions. As the eigenvalues remain purely imaginary, only the magnitude of the imaginary part is indicated. From α = 0 to α = 0.5, the error for low frequencies decreases, while it increases for high frequencies. All frequencies are under-estimated with respect to the exact values. For α = 0.5, every two eigenvalue pairs coincide at values that approximate the original spectrum best. This centered parameter choice is certainly more appropriate for a diffusive system. We will therefore, in a next step, analyze the spectral approximation of the Dirac structure for the parabolic case of constitutive equations and boundary conditions from the heat equation.

Conclusions

We introduced the weak form of the Stokes-Dirac structure with boundary inputs different causality. This Stokes-Dirac structure is the underlying geometric structure to represent power continuity in a port-Hamiltonian distributed-parameter system. At the example of a system of two conservation laws with canonical interdomain coupling, we described the mixed Galerkin discretization of the Stokes-Dirac structure in a general way. To obtain finite-dimensional approximate models in port-Hamiltonian form with the prescribed boundary inputs -as basis for the interconnection of multi-physics models, control design and simulation -we proposed power-preserving projections. These maps on the space of discrete bond variables (the mixed Galerkin degrees of freedom associated with flows and efforts) allow to define non-degenerate duality pairings on a new bond spaces, leading to finite-dimensional approximate Dirac structures. The latter admit several representations, one of them the input-output-representation. Port-Hamiltonian state space models are obtained, if in addition, the energy and the constitutive equations are approximated consistently. At the example of Whitney finite elements we demonstrated the discretization procedure and gave interpretations of the resulting discretization schemes.

The proposed method is, to the best of our knowledge, the first method which allows for a structurepreserving discretization of port-Hamiltonian distributed-parameter systems in more than one spatial dimension with a systematic treatment of different boundary inputs and the possibility to tune the discretized models between centered schemes and upwinding. The proposed family of approximation Dirac structures avoids a direct feedthrough in the state space model and the over-estimation of higher frequencies in the approximate spectrum, which is the case for the method presented in [START_REF] Golo | Hamiltonian discretization of boundary control systems[END_REF], where the efforts instead of the flow degrees of freedom are projected. The weak form of the Stokes-Dirac structure is the key feature that allows to include additional effects such as dissipation or diffusion or, more generally, to tackle the discretization of port-Hamiltonian systems with general and higher order interconnection operators and distributed inputs.

Remark 1 .

 1 Note that a Stokes-Dirac structure does not lead necessarily to the PH representation of a hyperbolic system of conservation laws as in the two examples above. With the corresponding definition of flow and effort variables and closure equations, parabolic phenomena based on irreversible conservation laws of thermodynamics (heat conduction) can be modeled based on the same canonical Stokes-Dirac structure. See the 3D heat conduction example in[START_REF] Duindam | Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach[END_REF] (Section 4.2.2, p. 233) or the relation[START_REF] Zwart | Building systems from simple hyperbolic ones[END_REF].

Proposition 1 .

 1 Given the n-dimensional open and connected domain Ω with Lipschitz boundary ∂Ω. Consider ∂Ω completely covered by subsets Γ i ⊂ ∂Ω, i = 1, . . . , n Γ , and Γj ⊂ ∂Ω, j = 1, . . . , nΓ with orientation according to ∂Ω. Let n Γ i=1 Γ i ∪ nΓ j=1 Γj = ∂Ω and the intersections Γ i ∩ Γj are either empty or have dimension n -1. Define the boundary flow and effort forms

2 .

 2 Ωlk denotes the common part of the boundary of Ω k and Ωl , respectively, where the minus sign underscores the inverse orientation by construction. For an illustration of the case n = 2, see Fig. 1. Define on each common boundary ∂Ω kl = -∂ Ωlk the interconnection conditions f kl k = e p | ∂Ω kl = êlk l and e kl k = (-1) p e q | ∂Ω kl = f lk l . Then, the terms e kl k |f kl k ∂Ω kl and f lk l |ê lk l ∂ Ωlk in the overall power continuity equation cancel each other out due to the reverse integration direction. The interconnection is hence power-preserving, and the composition of the separate Stokes-Dirac structures is a Stokes-Dirac structure due to their compositionality property.

Figure 1 :

 1 Figure 1: Sketch of a 3-dimensional domain Ω with subdomains Ω 1 , Ω 2 , Ω 3 and Ω1 and a subdivision of the boundary into Γ 1 , Γ 2 , Γ 3 and Γ1

Definition 4 .

 4 The vectors of discrete boundary port variables 25 e b,µ , f b,µ 0 ∈ R M µ b and êb,ν , f b,ν

  ) the resulting state space model ẋ = J∇H d (x) + Bu y = B T ∇H d (x) + Du (93) has PH form and the discrete energy satisfies the balance equation Ḣd = -y T u,

Figure 2 :

 2 Figure 2: Left: Numbering of nodes, edges and faces on a N × M rectangular simplicial mesh. K = (M + 1)N , L = K + M (N + 1). Right: 3 × 2 mesh.

Figure 3 :

 3 Figure 3: Illustration of an element of the row elements of P f p , associated to a discrete flow variable f p i . Defined as a weighted sum of original discrete flows on the neighboring 2-simplices, it can be considered also topologically dual to ẽp i .

Figure 4 :

 4 Figure 4: Weights on directed edges to define the elements of f q which are dual to ẽq . (a), (b) and (c) show the construction based on discrete efforts in horizontal, diagonal and vertical direction (red). The blue edges and weights display the contributions to the discrete constitutive equation, which is a flux across the effort edge. δ I/II = 1 8 + 1 4 (α I/IIβ I/II ), ε I/II = 1 8 -1 4 (α I/IIβ I/II ).

Figure 5 :

 5 Figure 5: Illustration of the elements of the output vectors f b = S p e p and f b = Ŝq e q as a weighted sum of node/edge efforts in the neighborhood of the input edge/node.

1 b

 1 ∈ R M b ×Nq is the matrix composed of unit row vectors associated to boundary edges, and the inputs to the simulation model êb (t) = u(t), e b (t) = 0. (102)

Figure 6 :

 6 Figure 6: Propagation of a wave due to point-wise boundary excitation under different parametrizations of the method. Snapshots at t sim = 16.

Table 1 :

 1 Dimensions of discrete flow and efforts spaces on the rectangular N × M simplicial grid.

	Vector(s) Dimension 2N M 3N M +N +M (N +1)(M +1) 2(N +M ) f p f q , e q e p e p b , e q b Symbol(s) Np Nq = Mq Mp M b p , M b q

Table 2 :

 2 Sizes and ranks of the matrices for the discrete variational problem and the direct discret model, respectively. N, M > 2.

	Matrix Size Rank	Mp Mp × Np Mp × Nq Np × Nq Kp + Lp dp Mp -2 Mp -2 Np	Mq Nq × Nq 2(Mp -2)	Kq + Lq Nq × Mp Nq × Mp dq Mp -1 Mp -1	Lp = L T q Mp × Nq 2(M + N ) -1

Table 3 :

 3 Parameter sets for the simulation.

	αI Parameters 1 1/2 Parameters 2 1/2 Parameters 3 2/3 Parameters 4 7/8 1/32 1/32 7/8 αII βI βII 1/3 1/3 1/2 1/4 1/4 1/2 1/6 1/6 2/3

Table 4 :

 4 Values of the imaginary parts for the discretized wave equation with different projection parameters, compared to the exact values. Grid size N = 20. .8451 26.6130 24.4593 22.6179 17.4897 15.0650 12.9890 20 61.2611 39.8826 31.9591 24.9858 21.8185 20.4059 19.9383

	k 1 2 3 4 10.9956 10.5993 10.5115 10.4509 10.4193 Exact α = 0 α = 0.1 α = 0.2 α = 0.3 α = 0.4 α = 0.5 1.5708 1.5321 1.5392 1.5464 1.5539 1.5614 1.5692 4.7124 4.5873 4.5985 4.6122 4.6284 4.6473 1.5692 7.8540 7.6156 7.6013 7.5974 7.6045 7.6231 4.6689 8.4371 4.6689 5 14.1372 13.5207 13.2958 13.1256 13.0148 9.6285 7.6537 10 29

The skew-symmetric wedge product ∧ : Λ k (Ω) × Λ l (Ω) → Λ k+l (Ω) pairs two differential forms, adding their degrees.

See[START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF] for several examples

For brevity, all arguments are omitted

In 2D we have * dx = dy, * dy = -dx.

If not,[START_REF] Bossavit | Generating whitney forms of polynomial degree one and higher[END_REF] still represents a PH system, as the rotational term does not contribute to the energy balance. It can be associated to a canonical Stokes-Dirac structure, with a different definition of the dynamic equation for the 1-form u .

As a reference to bond graph modeling of dynamical systems[START_REF] Henry | Analysis and design of engineering systems[END_REF], see also[START_REF] Duindam | Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach[END_REF], Chapter 1)

or another field

For brevity, the domain of the differential forms is written as an index, Λ p Ω = Λ p (Ω).

This is definitely true for finite elements, which are the standard Galerkin approximation in more than one spatial dimension. Also for the 1D pseudo-spectral collocation method, see[START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF].

Discrete boundary variables have index b, in contrast to index ∂ for the original distributed quantities.

This denomination refers to the trace theorem that clarifies the functional space of H 1 functions restricted to (a subset of) the boundary.

This is the case of a relaxed image representation.

If a flow variable is defined differently, as in the case of the 2D SWE with the additional rotation term, this has to be accounted for also in the discrete equation.

The weights can be different from node to node. For simplicity, we fix the same pairs (α I , β I ) and (α II , β II ) on the whole mesh.

Precisely, the average value of the coefficient function of the 2-form p.

Index raising of the 1-form q.

Note that the plots in Fig.6represent the discrete, projected efforts ẽp i in the nodes of the mesh. An appropriate localization of the values of p would be the weighted barycenters of the areas according to Fig.3. This would advance the plotted wavefront in the propagation direction.

Current and future work concerns the application of the method to the port-Hamiltonian representation of the heat equation, which shares the same Stokes-Dirac structure, the analysis of system-theoretic properties of the discretized models in view of control design and the interface to existing finite-element and simulation software. Moreover, we intend to include the discretization of the nonlinear constitutive relations for the 2D shallow water equations in our open models and clarify the links with recent work on geometric mixed finite elements like [19], [20], where in-and outputs are not explicitly taken into account, and upwinding in differential forms as presented in [18].