{'fraction_non_alphanumeric': 0.036846967053269016, 'fraction_numerical': 0.012888374070761426, 'mean_word_length': 4.450927109974424, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dual field structure-preserving discretization of port-Hamiltonian systems using finite element exterior calculus

In this paper we propose a novel approach to discretize linear port-Hamiltonian systems while preserving the underlying structure. We present a finite element exterior calculus formulation that is able to mimetically represent conservation laws and cope with mixed open boundary conditions using a single computational mesh. The possibility of including open boundary conditions allows for modular composition of complex multi-physical systems whereas the exterior calculus formulation provides a coordinate-free treatment. Our approach relies on a dual-field representation of the physical system that is redundant at the continuous level but eliminates the need of mimicking the Hodge star operator at the discrete level. By considering the Stokes-Dirac structure representing the system together with its adjoint, which embeds the metric information directly in the codifferential, the need for an explicit discrete Hodge star is avoided altogether. By imposing the boundary conditions in a strong manner, the power balance characterizing the Stokes-Dirac structure is then retrieved at the discrete level via symplectic Runge-Kutta integrators based on Gauss-Legendre collocation points. Numerical experiments validate the convergence of the method and the conservation properties in terms of energy balance both for the wave and Maxwell equations in a three dimensional domain. For the latter example, the magnetic and electric fields preserve their divergence free nature at the discrete level.


Introduction
The field of computational engineering seeks to advance the reliability and accuracy of simulations by means of faster and/or higher fidelity algorithms. Complex engineering systems arise from the interconnection of many simple components. For this reason, it is important that computational models be also constructed via a modular approach, where each subsystem is modelled separately and then assembled to the rest (as in dynamic substructuring). Furthermore, even more important is the fact that such models can be reusable and are to a big extent independent of the super-model they will be part of. This idea is at the core of the port-Hamiltonian (pH) framework, since it describes the interaction between systems and the surrounding environment by means of ports (cf. [1] for the first paper discussing pH finite dimensional systems and their connection with bond-graph modelling). In the infinite dimensional case the mathematical literature has been mostly focused on Hamiltonian systems with no interaction with the environment (see [2] for a reference on the subject). However in the seminal work [3], distributed pH systems were first introduced to model physical systems of conservation laws (which happen to be described by hyperbolic PDEs due to the physical finite propagation time). In order to account for the non-trivial boundary interaction, the authors introduced the concept of Stokes-Dirac structure.
This geometrical structure represents an extension of Dirac manifolds [4] to the infinite dimensional case. The Stokes theorem allows defining appropriate boundary variables to account for the power exchange through the boundary of the spatial domain. This geometrical structure is acausal, in the sense that it ignores the actual boundary conditions of the problem at hand, and simply captures all admissible boundary flows.
Distributed pH systems bring together different areas of mathematics, physics and engineering. They were first formulated using the language of differential geometry but their mathematical as well as system theoretic properties have also being studied in a functional analytic setting (cf. [5,6,7,8] for works on mathematical wellposedness and [9,10] for studies on stability and stabilization). Since they are closed under interconnection (this feature stems from the properties of the Dirac structure [11]), pH systems have the potential to tackle complex multiphysical engineering applications. So far they were employed to model fluid-structure coupled phenomena [12], reactive flows [13], Euler and Navier-Stokes equations [14,15,16], thin mechanical and thermomecanical structures [17,18]. The interested reader may consult [19] for a comprehensive review on distributed pH systems.
To exploit the full potential of this novel modelling framework, discretization schemes must retain the many properties of pH systems. In particular, numerical algorithms should capture the underlying conservation laws and provide a systematic way to handle boundary conditions, so that the resulting discrete pH system defines a (modular) Dirac structure. The first contribution in this sense dates back to [20]. Therein, the Stokes-Dirac structure is discrete in strong form exploiting mixed finite element differential forms and compatibility conditions need to be fulfilled 1 . A reduction method based on symplectic Gauss-Legendre collocation points is detailed in [22]. This methodology suffers a major drawback, as it cannot handle higher dimensional spatial domains. A staggered finite difference strategy is proposed in [23] but the scheme is specialized for the wave equation only. The discrete exterior calculus framework [24], based on algebraic topology, can also be used to formulate pH systems in a purely discrete manner, as shown in [25]. However, this framework does not explicitly rely on interpolating basis function but only on mimetic discrete operators that reproduce the behaviour of their continuous counterparts. This makes necessary the employment of dual topological meshes, based on the Delaunay-Voronoi duality, to construct an isomorphic discrete Hodge star. This limits the applicability of the method, as dual meshes require ghost points that do not lie inside the physical domain, thus making it difficult to interconnect systems along a physical boundary. In order to avoid the need for dual meshes, a Galerkin formulation based on Whitney forms is detailed in [26]. In order to construct a minimal bond space 2 , projector matrices parametrized by tunable parameters need to be introduced. The best choice of these parameters depends strongly on the application at hand and this complicates the employment of the method. When only one computational mesh is used, it is not possible to construct a discrete isomorphic Hodge operator [28,29]. Mixed finite elements do not rely on dual meshes to construct an isomorphic discrete Hodge, but rather on a weak formulation of the codifferential operator via the integration by parts formula. This was formalized by the mimetic discretization community [30,29,21] and, in the context of port-Hamiltonian systems, by the Partitioned Finite Element method [31]. In particular, in [32] convergence of mixed finite element discretization is guaranteed for a large class of finite elements, and not only restricted to the families that satisfy the de Rham cohomology.
From the past literature, it may look like new tools and implementations are needed to be able to capture the particular features of pH systems. In this paper we argue the opposite. The pH modelling paradigm represent a unifying framework to describe physical problems. As such, it should be discretized via a unifying discretization method: the Finite Element Exterior calculus theory, initiated by the Japanese papers of Alain Bossavit and fully developed by Arnold, Falk and Whinter [21], provides all the tools needed to accomplish this task. Using finite element differential forms, a discrete counterpart of the Stokes-Dirac structure can be easily obtained. However, the discretization of constitutive equations remains problematic as those require a discrete Hodge star. Recently, a dual field formulation has been introduced to discretize the Navier-Stokes equations [33]. This methodology seeks variables in dual polynomial finite element spaces via a mixed finite element formulation and is capable of conserving mass, helicity and energy. A staggered time discretization, in which the vorticity information is exchanged between the two systems, allows treating the non linear rotational term in a linear fashion and makes this scheme computationally competitive with state of the art methods. This method has a clear connection with Hamiltonian dynamical systems, as discussed by the authors in the introduction of the paper. Furthermore, it relies on the Finite Element Exterior calculus framework as the finite elements used therein (the so-called mimetic quadrilateral/hexahedral finite elements [34]) form a discrete de Rham complex.
In this contribution, we show how the dual field formulation can be used for the systematic discretization of linear pH systems (but the framework is extendable to the state modulated case where one has a Poisson manifold, as already done in [33]). To show the generality of the approach, we detail the calculations in exterior calculus. The dual field formulation is based on the construction of an adjoint Stokes-Dirac structure, where the flows and efforts are the Hodge dual of the initial Stokes-Dirac structure. The introduction of an adjoint system leads to a redundant representation of the dynamics, as the state variables are doubled. Nevertheless, this redundancy removes the need for a discrete Hodge star, as one can use dual variables to expresses the energy as a duality product of forms, rather than the classical inner product. The Stokes-Dirac structure can be immediately discretized in a mimetic fashion, whereas the adjoint Stokes-Dirac structure requires the integration by parts formula to obtain a discrete representation of the codifferential. The constitutive equations are embedded in the dynamics by expressing the state variables in terms of the effort ones. In other words, by using a redundant dual system, the state variables will be the original ones, plus the Hodge star of the former. In this way, it is possible to calculate the needed efforts corresponding to the differential of the energy function, without any need of the Hodge operation. The boundary conditions are here imposed strongly. This choice leads to two uncoupled mixed finite element discretization of the dynamical pH systems, each containing two unknowns (as in the Partitioned Finite Element method). In this way, after these two systems are integrated on the same time grid via symplectic collocation Runge-Kutta methods (for simplicity the implicit midpoint scheme will be used here), an exact discrete power balance is retrieved. This implies that the discretization scheme gives rise to a Dirac structure and can be used for interconnecting different physical systems in a structured fashion. We detail the discrete operations (like the inner and duality products, the exterior derivative and the trace) using the trimmed polynomial spaces. The lowest order case coincide with the Whitney forms [35]. The Whitney forms are particularly meaningful, since they make explicit a number of concepts like metrical versus topological properties and operations. Numerical tests, assessing the conservation properties of the scheme and the convergence rate of the different variables, are performed on the wave and Maxwell equations in a three dimensional Euclidean domain 3 . The numerical tests show that each of two mixed discretizations converge to the exact solution under mesh refinement, without exhibiting synchronisation problems. The well established open-source library Firedrake [37] is employed for the implementation (but one can also use the FEniCS [38] library).
The paper is organised as follows: in Sec. 2 some preliminaries on needed concepts and notation will be presented. In Sec. 3 we introduce the concept and representation of the classical and adjoint Stokes-Dirac structure and the associated pH systems which is instrumental in the presented solution. The weak formulation, upon which the discretization scheme is based, is presented in Sec. 4. The spatial discretization is illustrated in Sec. 5 by means of the trimmed polynomial spaces. We detail the discrete version of the dual field Stokes-Dirac structure as well as the associated dynamical system. In Sec. 6, the time discretization is explained using the implicit midpoint method. Sec. 7 collects the numerical tests and shows the numerical properties of the proposed methodology.

Preliminaries
In this section we will introduce some preliminary concepts and their notation which are needed for the treatment ahead. The presented material is inspired by [29,36].

Smooth differential forms
Let M ⊂ R n be a differentiable Riemannian manifold with boundary ∂M and metric g. The space of smooth differential forms on M (i.e. the space of smooth sections of the kth exterior power of the cotangent bundle T * M ) is denoted by Ω k (M ). The wedge product ∧ : Ω k (M ) × Ω l (M ) → Ω k+l (M ), k, l ≥ 0, k + l ≤ n is the skew-symmetric exterior product of differential forms One fundamental operator acting on smooth differential forms is the exterior derivative d :

that satisfies the following axiomatic properties
• df for f ∈ Ω 0 (M ) is the differential of f ; Occasionally for the sake of additional clarity, the exterior derivative acting on a k-form will be denoted by d k . Consider a manifold M of dimension n, a manifold S of dimension m, and a smooth mapping between them, Φ : S → M . The pullback operator, Φ * , is a mapping Φ * : Ω k (M ) → Ω k (S), that maps k-forms in M to k-forms in S, with k ≤ m, and k ≤ n, naturally. An important case of pull back is the trace operator.
Definition 1 (Trace operator). The trace operator is defined to be the pull back of the inclusion map ι : One fundamental result in exterior calculus is the Stokes theorem, that relates derivation, integration and the trace operator with M an n−dimensional manifold. Differential forms possess a natural duality product.
Definition 2 (Duality product). Given a smooth manifold M of dimension n, the duality product is denoted by The duality product is also defined on the boundary ∂M (whose orientation is inherited from the one of the manifold) where k = 0, . . . , n − 1. Combing together the Leibniz rule and the Stokes theorem, one has the integration by parts formula 2.2. L 2 theory of differential forms Given a coordinate chart ξ i : M → R i = 1, . . . , n, we can represent locally a point p in the manifold M with the tuple ξ(p) := (ξ 1 (p), . . . , ξ n (p)). Then, the local representation of a form where the multi-index I = i 1 , . . . , i k , 1 ≤ i 1 ≤ · · · ≤ i k ≤ n has been introduced. The space of k-forms can be equipped with a pointwise inner product, inherited form the metric structure of the Riemannian manifold (that establishes an inner product of vectors) and the duality between vectors and differential forms.
Definition 3 (Pointwise inner product of forms). Let α k , β k ∈ Ω k (M ) be two forms with local representation as in (8), the point-wise inner product is given by where g k,l = (dξ k , dξ l ) are the components of the inverse metric tensor.
Once an orientation is given to the Riemannian manifold, the Hodge star operator can be properly defined. The Hodge maps inner-oriented (or true) forms, measuring intensities, to outeroriented (or pseudo) forms [39,40], measuring quantities, and vice versa. This distinction is of fundamental importance and allows defining quantities that are orientation independent (like mass, energy, etc.). In this paper outer-oriented forms are denoted by means of a hat, i.e. α k ∈ Ω k (M ) where Ω k (M ) is the space of outer oriented (or pseudo) forms.
Definition 4 (Hodge-operator). The Hodge-operator, defined for an n-dimensional Riemannian manifold M , is the operator : Ω k (M ) → Ω n−k (M ) such that where the inner product is defined in Def. 3. The standard volume form in local coordinates is given by [40, Page 362] where g ij are the components of the metric tensor in the chosen chart ξ and 1...n is the Levi-Civita symbol. Notice the often forgotten importance of the Levi-Civita symbol, that indicates that the volume form is indeed a pseudo form vol = 1, as the total volume of a portion of space cannot be negative.
We now introduce the L 2 inner product of forms.
Definition 5 (L 2 inner product). Given a smooth manifold M of dimension n, the L 2 inner product is defined by As in vector calculus, the L 2 Hilbert space of differential forms is defined by completion.
Definition 6 (L 2 space of differential forms). The L 2 Ω k (M ) space is the completion of the space of smooth forms Ω k (M ) in the norm induced by the L 2 inner product.
The L 2 inner product is fundamental for the construction of weak formulations. Weak formulations allows enlarging the solution space from the space of smooth forms, to the less regular Sobolev spaces, that represent the domain of the exterior derivative.
Definition 7 (Sobolev spaces of differential forms). The space of square integrable forms whose derivative is square integrable is denoted by For k = n, the space reduces to an L 2 space, HΩ n (M ) = L 2 Ω n (M ).
These spaces, connected by the operator d, form the de Rham domain complex, as illustrated in the following diagram. The adjoint 4 operator to the exterior derivative with respect to the L 2 inner product is the codifferential.
The definition is such that the codifferential corresponds to the formal adjoint of the exterior derivative 5 Sobolev spaces for the codifferential are then defined as the space of forms whose codifferential is square integrable. In this work the codifferential will be always treated in a weak manner, using the integration by parts in Eq. (11). This choice leads to a weak formulation in which all variable live in the de-Rham complex of Fig. 1.

The canonical and adjoint Stokes-Dirac structure
In this section the strong formulation of port-Hamiltonian systems and the important role of the underlying geometrical structure are recalled. The corresponding adjoint system is constructed. For sake of concreteness, first the introductory example of the scalar wave equation with Dirichlet boundary control is presented. This problem is a simple example of port-Hamiltonian system.

An introductory example: the scalar wave equation
The propagation of acoustic waves in M ⊂ R 3 is described by the following hyperbolic partial differential equation, that determines the time-dependent (pseudo) 3-form field η 3 (t) : together with time-varying Dirichlet boundary condition To highlight the Hamiltonian structure of the wave equation, consider the variables Equation (12), together with the boundary condition (13), can now be recast into a first order system Here u 0 corresponds to a control input applied to the boundary. Indeed port-Hamiltonian systems are boundary controlled systems and the boundary conditions coincide with inputs that describe interactions with the external environment.
The choice of variables v 3 , σ 1 leads to a clear separation between topological and metrical properties. The topological properties are incorporated into the (skew-dual) operator J, whereas the metric properties are evident in the appearance of the Hodge in the right hand side of Eq. (12). The topological properties relate to conservation laws, where the metrical properties are related to the constitutive laws of the physical system under consideration.
The power that flows through the domain boundary is obtained by taking the duality product with v 3 and σ 1 and using the Leibniz rule and Stokes theorem where variable y 2 corresponds to the power conjugated output to the input, i.e. the power collocated output y 2 := − tr σ 1 .
It can be noticed that this variable corresponds to the Neumann boundary condition [36]. Furthermore, the output is dual to the input on the boundary ∂M . Notice that each duality product in (16) involves pairing an inner-oriented form with an outer-oriented one.
System (15) is an example of a port-Hamiltonian system. The underlying geometrical structure is the Stokes-Dirac structure, an infinite dimensional generalization of Dirac manifolds introduced by Courant [4]. In what follows we introduce the general port-Hamiltonian formulation for a system of two conservation laws.

Linear Stokes-Dirac structures
Given a smooth manifold M of dimension n, the flow variables f p 1 ∈ Ω p (M ), f q 2 ∈ Ω q (M ) and the effort variables e n−p 1 ∈ Ω n−p (M ), e n−q 2 ∈ Ω n−q (M ), with 6 p + q = n + 1, consider the Stokes-Dirac structure where r = pq + 1 for mathematical and physical reasons [3]. The operator J, representing what is called in physical system theory the junction structure, is called the interconnection operator and is a formally skew-dual operator. This property is at the core of the balance equation where ∂M denotes the boundary of manifold M . This balance equation arises from the application of the Leibniz rule and the Stokes theorem and states the overall conservation of energy. This is the embodiment of Tellegen's theorem generalised to the distributed parameters systems case [41].
Remark 1. The Stokes-Dirac structure is a geometrical concept that does not involve time. The exterior derivative is associated to the spatial manifold only. A flow variable may be associated to a dynamics (i.e. f k = −∂ t α k as it will be shown in the next section) leading to a conservation law, or to a resistive relation that introduces dissipation. The hyperbolic case arises when both flows have an associated dynamics, whereas the parabolic case is obtained when one flow is associated to a dynamics and the other one is purely algebraic and associated to resistive relation. In this paper the focus will be on the hyperbolic case only.

Port-Hamiltonian systems
Port-Hamiltonian systems encoding conservation laws are associated to the geometrical Stokes-Dirac structure. To establish this connection, consider the distributed state variables α p (t) : → Ω q (M ) and the Hamiltonian functional with Hamiltonian density n-form H. A fundamental notion is the variational derivative of the Hamiltonian functional [2,3].
Definition 9 (Variational derivative). The variational derivatives of the Hamiltonian δ α H n−p ∈ Ω n−p (M ), δ β H n−q ∈ Ω n−q (M ) are defined implicitly by Considering then a trajectory ( α p (t), β q (t)) parameterised by time in the manifold of state fields, the variational derivative is defined so that the rate of the Hamiltonian (20) readṡ With an abuse of notation, we will write α p (ξ, t) = α p (t), β q (ξ, t) = β q (t) to make explicit the value of the field at a certain point ξ ∈ M rather than the field as a section. Consider then the following equations which can be seen to represent a system of two conservation laws with canonical inter-domain coupling [3] with initial conditions The flows and efforts of the associated Stokes-Dirac structure (18) are then defined by Boundary conditions. Mixed boundary conditions will be considered in this work. We denote with Γ 1 and Γ 2 two open subsets of the boundary that verify Γ 1 ∪ Γ 2 = ∂M and Γ 1 ∩ Γ 2 = ∅. Each boundary subpartition is associated with one boundary condition. In particular, the values of δ α H n−p = e n−p 1 (resp. δ β H n−q = e n−q 2 ) are imposed on Γ 1 (resp. Γ 2 ). Since in this paper we consider boundary controlled systems, the boundary conditions are assigned by means of the inputs tr e n−p The outputs are defined in such a way that they are power conjugated to the inputs Given the balance equation (19), the rate of change of Hamiltonian (power flow) readṡ where the last equality descends from the additive property of the integration and the definition of the boundary variables given in Eq. (18). Equation (27) states that the change of energy within the spatial domain M is equal to the supplied power through its boundary ∂M . Notice that since p + q = n + 1, the boundary variables are dual forms on the boundary (n − p) + (n − p) = 2n − (n + 1) = n − 1.
Constitutive equations. Since in this work the focus is on the case in which the closure equations are linear, we make the additional assumption of a quadratic Hamiltonian.
Assumption 1 (Quadratic Hamiltonian). The Hamiltonian density in assumed to be quadratic in the state variables where A p : L 2 Ω p (M ) → L 2 Ω p (M ) and B q : L 2 Ω q (M ) → L 2 Ω q (M ) are bounded from above and below, symmetric, positive definite operators with respect to the standard inner products in L 2 . These tensors are related to physical properties of space like electric permittivity, mass density or Young modulus.
Because of this assumption, the variational derivatives of the Hamiltonian, computed according to Definition 9, can be seen to be

Adjoint Stokes-Dirac structure
The port-Hamiltonian system (22) and its associated Stokes-Dirac structure (18) were defined on the state space denoted by X = Ω p (M ) × Ω q (M ). The key ingredient that will be used for our proposed dual-field discretization approach is the Adjoint Stokes-Dirac structure that we will introduce in this section. The adjoint structure is defined on the Hodge dual state space It is worth noting that the choice of Φ is not unique. While in this work we choose (30) for numerical reasons, in Sec. 8 we will discuss an alternative in which the material tensors will be included in the Hodge. The pushforward associated to the diffeomorphism converts elements in the tangent space of the state space the adjoint Dirac structure of (18) is expressed by means of co-differential map d * (10) as where the following notation has been used a 0 = r + p(n − p) + q(n − q) + n(q + 1) + 1, a 1 = n(p + 1) + 1.
The following corollary is an immediate consequence of Prop. 1.

Corollary 1.
The adjoint Stokes-Dirac structure verifies the power balance Proof. The property follows from the fact that the pullback (32) preserves duality products. It can be verified by a direct computation that: Symmetry of the inner product, , Property of the the Hodge star, Analogous computations are performed for the second flow and effort variables f q 2 and e n−q 2 .
To simplify the expressions of the adjoint Dirac structure and relate the new indices a 0 and a 1 to the p, q, n, the following proposition is of value.
Proposition 2. The coefficients a 0 , 1 and a 1 , 1 + r + p(n − p) + q(n − q) have the same parity which means a 0 ≡ 1 mod 2, where mod 2 denotes the modulo of the division by 2.

Proof. Reported in Appendix A.
This results implies that: and this could be used to express (35) with the coefficients p, q, n and r, and will be of importance later. The adjoint Stokes-Dirac structure introduces the codifferential and its important role in discretization.

Associated adjoint port-Hamiltonian system
Given the adjoint Stokes-Dirac structure (35) and the pH system (22), the adjoint pH system is obtained by considering the isomorphism and considering the new Hamiltonian as a function of the state variables of the adjoint system The energy density H corresponds to H under the isomorphism Φ. The variational derivative, being defined by duality product (cf. Def. 9), undergoes the same pullback as the efforts, so that The port-Hamiltonian adjoint system is obtained substituting into the adjoint Stokes-Dirac structure (35) the dual effort variables (41) and the dual flow variables Given Assumption 1, the adjoint energy is a quadratic function in the adjoint state variables where the operators A n−p , respectively B n−q are the dual of A p , respectively B q defined implicitly by are bounded from above and below, symmetric and positive definite with respect to the L 2 inner product. Furthermore they verify the commutativity property A p = A n−p , B q = B n−p , thus are explicitly expressed by The variational derivative of the adjoint energy then reads For what concerns the boundary conditions, it is important to notice that the adjoint system carries the same boundary conditions (25) as the canonical one, but expressed by means of the adjoint efforts tr e p and analogously for the outputs.
Remark 4. By comparing the adjoint energy density (44) to the original one in (28) it can be seen that the adjoint port-Hamiltonian system (43) describes the dynamics governing the true-form representation of the α state and the pseudo-form representation of the β state. On the other hand, the canonical port-Hamiltonian system (22) describes the dynamics of α as a pseudo-form and β as a true-form. Combining the two representations is the key idea behind the dual-field formulation which we discuss next.

Weak dual field formulation
We now present the weak dual field formulation of the port-Hamiltonian model, inspired by the work of [33], but now including general mixed boundary conditions. In this formulation, we combine together the canonical and adjoint port-Hamiltonian systems described in the previous section. This combination will eliminate the necessity of using an explicit Hodge star operator, which will be accounted for using the integration by parts formula in the weak adjoint system (11). As a consequence, two decoupled port-Hamiltonian systems are obtained: one based on outer-oriented forms, and a second, based on inner-oriented forms. In the following we will refer to the former as the primal system and at the latter as the dual system.

Motivation of the approach
The novelty of this approach is that the Hamiltonian can now be rewritten as a function of the original state variables and their Hodge duals: with associated density The two expressions for the Hamiltonian density are equivalent at the continuous level but not when one seeks a discrete representation of the variational derivative. In the second expression, the Hodge star operators have been removed and effectively replaced by a dual-field representation of the variables. This is important because the discretization can be performed without relying on a discrete Hodge star. The effort variables are now computed using the latter expression of the total Hamiltonian The constitutive equations are then incorporated directly in the dynamics by reducing the state variables in terms of the efforts (in the literature this is called co-energy formulation [17]). This reduces the number of equations to be solved and formulate the problem in terms of the most important variables, as the efforts are subjected to the boundary conditions. To express the dynamics in term of the effort variables the material operators A p , B q , A n−p , B n−q are inverted (this is possible since these operator are positive definite and bounded from above and below for physical reasons) where the tensors These tensors are commonly introduced in mixed finite element (for example in linear elasticity they represent the density and the compliance [43], whereas in electromagnetism they represent the electric permittivity and the magnetic permeability [44]).

Dual field Stokes-Dirac structure
Now we will present the dual field Stokes-Dirac structure which will be constructed by combining the Stokes-Dirac structure (18) and its adjoint (35). In order to obtain the weak form, we take the inner product of (18) and (35) with test differential forms of appropriate degrees and use the integration by parts formula (11) for the adjoint system. This gives rise to two uncoupled systems. One formulation contains outer-oriented forms only, except for the boundary term arising from the integration by parts. We will refer to this formulation as the primal formulation. The second formulation contains inner-oriented forms only, except once again for the boundary term. We will refer to this formulation as the dual formulation.
Primal Stokes-Dirac structure. The primal Stokes-Dirac structure contains outer oriented forms. Considering that p + q = n + 1, it holds n − q = p − 1 and n − p = q − 1. This shows that the primal formulation corresponds to a first mixed discretization in which the variables are related by d p−1 (52) where in the integration by parts, the fact that H −1/2 Ω q−1 (∂M ) f q−1 ∂ = tr e q−1 1 = tr e p 1 has been used. The fractional Sobolev spaces H −1/2 Ω q−1 (∂M ), H −1/2 Ω p−1 (∂M ) arise from the fact that the trace operator on H 1 Ω k (M ) (the space of k-forms whose coefficients are in H 1 (M )) extends to bounded linear operator tr : The dual Stokes-Dirac structure. The dual Stokes-Dirac structure contains inner oriented forms. This formulation corresponds to a second mixed discretization in which the variables are related by d q−1 where it has been used a 0 ≡ 1, mod 2 and tr e q 2 = (−1) q(n−q) tr e p−1

Boundary conditions enforcement
The primal and dual Stokes-Dirac structures have an associated primal and dual port-Hamiltonian system respectively. The Stokes-Dirac structures do not carry boundary conditions information but the associated port-Hamiltonian systems do. A number of different techniques can be used to enforce boundary conditions (cf. [45] and references therein). We opt for a strong imposition of the essential boundary condition (i.e. the one that does not appear in the integration by parts) by direct assignment of the boundary degrees of freedom in each resulting port-Hamiltonian system. In this way two uncoupled port-Hamiltonian systems are obtained. Furthermore, this strategy will ensure that the power flow is correctly represented (up to the approximation error introduced by the polynomial interpolation) as it will be proven later.

Strong primal and dual port-Hamiltonian system
In this section, the primal-dual port-Hamiltonian system will be constructed. For sake of clarity the systems are first presented in strong form and then the associated weak form is presented.
Strong primal port-Hamiltonian system. Using the first line of (22) and the second line of (43), together with the expression of the state variables in terms of the effort in Eq. (51) and Pr. 2, the outer oriented primal system is obtained: where e p 1 ∈ Ω p (M ) and e p−1
Strong dual port-Hamiltonian system. Using the first line of (43) and the second line of (22), together with the expression of the state variables in terms of the effort in Eq. (51), the inner oriented dual system is obtained: where e q−1

Weak primal and dual port-Hamiltonian system
The weak formulation for the primal and dual system is here presented. The material tensors are supposed to be positive and bounded above and below operators from L 2 Ω k (M ) to L 2 Ω k (M ) (analogously to what is done in [21,Sec. 7.3]) and are therefore allowed to be discontinuous. For this reason, a weak formulation of the dynamics is necessary. The codifferential is here interpreted weakly. This will allow incorporating the natural boundary condition as a system input.
The weak primal port-Hamiltonian system. The weak formulation for the primal system reads: The weak dual port-Hamiltonian system. For the dual system, the following weak formulation is obtained: find e q−1 (59)

Space discretization
In this section, the space discretization of the dual-field port-Hamiltonian model is performed using the trimmed polynomial family (whose construction is discussed in [21,36]). First, we start by introducing the discrete operators on the finite-element space which will mimic their continuous counterparts. The crucial point is that the employed finite element forms form a discrete de Rham complex, as stated in [21,Lemma 3.8].
For notational simplicity the trimmed polynomial space will be denoted by V k s,h = P − s Ω k (T h ), where k is the degree of the differential form, s the polynomial degree associated to the interpolating bases functions and h the mesh size. Furthermore, the notation [A] j i will be used to denoted the element of matrix A corresponding to the i-th row and the j-th column.

Mimetic operations
For a generic k-discrete form µ k h ∈ V k s,h , one has where N k,s is the number of degrees of freedom for V k s,h , the degree of freedom µ k i ∈ R is given in [21,Eq. 5.2], and φ k s,i : M → V k s,h ⊂ HΩ k (M ) is a section of V k s,h , corresponding to a finite element basis function. For an outer oriented form µ k h ∈ V k s,h the following notation is used for the finite element expansion Inner product. Given two discrete forms (ν k h , µ k h ) ∈ V k s,h , the inner product reads where ν k , µ k ∈ R N k,s are the vectors collecting the degrees of freedom µ k i , µ k i respectively and the mass matrix M k s ∈ R N k,s ×N k,s of order k (symmetric and positive definite) is computed as Duality product over the domain. Given two discrete forms variables λ n−k h ∈ V n−k s,h , µ k h ∈ V k s,h , their duality product over the domain is expressed as where the matrix L k s ∈ R N n−k,s ×N k,s is computed via [L k s ] j i = φ n−k s,i |φ k s,j M . Since the wedge product, used in the definition of dual product introduced in Eq. (5), is alternating, it could be seen that the following holds where the matrix L n−k Exterior derivative. The expression of the exterior derivative is here specialized for the inner and duality product. Given a form ν k+1 Remark 5. The duality product over the domain is not used by the numerical scheme. However, the scheme is such to preserve the integration by parts formula (7) that contains it.
Trace. For a discrete differential form of order k ≤ n − 1 the trace operator is also defined.
Since the basis forms are defined locally for each element of the mesh, when the trace of a form is considered the simplices that do not lie on the boundary can be discarded in the expansion. Denoting the boundary of the simplicial complex by where N ∂ k,s = n−1 j=k #∆ j (∂T h ) is the number of all j-simplices (with k ≤ j ≤ n − 1) along the boundary for a polynomial differential form of order k. The degrees of freedom along the boundary µ k l,∂ are associated to j−simplices σ l j,∂ , k ≤ j ≤ n − 1 lying on the boundary. The trace matrix simply collects them Duality product over the boundary. Given two discrete forms variables λ h ∈ V n−k−1 s,h , µ k h ∈ V k s,h , their boundary duality product is expressed as where the matrix L k s, Using the trace operator relation (68), the duality product can be expressed using the basis functions that lie on the boundary, i.e. ψ k s,l where the boundary matrix Ψ k s, ) . Using Eq. (68), the following matrix relation is obtained Discrete integration by parts. The integration by parts formula (7) is also valid when the differential forms and the domain are not smooth (cf. [36, Theorem 6.3]) where H 1 Ω k (M ) is the space of k-forms with coefficients in H 1 (M ). Conforming finite element V k s,h ⊂ HΩ k (M ) do not possess the H 1 Ω k (M ) regularity (except for the case k = 0). However, the integration by parts holds at the discrete level when conforming finite element spaces are used This is proven considering that the Stokes theorem holds for the each element of the mesh and it extends to the whole mesh as all terms arising on inter-cell boundaries will cancel due to the special continuity properties of discrete differential forms. A proof of this statement based on vector calculus is reported in Appendix A for the wave and Maxwell equations in 3D. A formal proof of the integration by parts formula in an exterior calculus setting is out of the scope of this paper. Since Eq. (74) is valid ∀µ h , ∀λ h , the algebraic form of the Stokes theorem is obtained where the alternating property of the duality product has been used for the first term.
Hodge star operator. As already discussed, in this paper a discrete Hodge star operator will never be employed thanks to the dual-field formulation of the system dynamics. The following discussion will highlight the benefits of not using a discrete Hodge star which is in fact a non trivial task [28]. Essentially, an isomorphic discrete Hodge star operator requires dual meshes. The explanation for the statement resides in the fact that finite element differential forms of dual type do not possess in general the same number of degrees of freedom. For example, in the Whitney forms case each k-order finite element has a dimension equal to the number of k-simplices. In general, the number of k-simplices will be different than the number of n − k simplices. Hence it is not possible to find a bijection relating a finite element k-form and its n − k Hodge dual. To retrieve a discrete isomorphism, dual meshes based on the Voronoi-Delaunay duality are normally used [24,28]. However, in some discrete formulations this approach result in mesh entities located outside the physical domain, making it cumbersome to implement boundary interconnections of different systems, which is a fundamental feature of port-Hamiltonian systems that we aim to preserve at the discrete level.
If one mesh is used, a non-isomorphic discrete Hodge star can be constructed relying on a weak formulation. Given µ k h ∈ V k s,h and its dual variable µ k h = µ n−k h ∈ V n−k s,h a discrete Hodge can be constructed using a weak formulation of the equation leading to the two following algebraic realization It can be noticed that this discrete Hodge corresponds to a projection between dual polynomial space and therefore it inevitably introduces an error. This fundamental point is the main motivation behind the dual field method. The integration of the primal and adjoint system allows reconstructing primal and dual variables without relying on a discrete Hodge star. Instead, it is embedded into the adjoint system by means of the codifferential operator.
The Whitney form case. Finally, we conclude this section by considering the special case of Whitney forms i.e. lowest-degree polynomials (s = 1). The special case of the Whitney forms allows for a clear separation of topological and metric operations: the resulting discrete Stokes-Dirac structure is completely topological, whereas the adjoint Stokes-Dirac systems embeds the metric in the discrete representation of the codifferential [29]. For a generic k-Whitney form µ h ∈ V k 1,h := W k h , one has where N k = #∆ k (T h ), ∀k = 0, . . . , n is the number of k-simplices in the mesh, µ i ∈ R corresponds to the cochain coefficient associated with the simplex σ i k and the function w k i (ξ) := φ k 1,i (ξ) is the k-th order Whitney function related to the k-simplex indexed by i and function of the point manifold ξ. The particularity of the Whitney form case is the exterior derivative, since only in this case the interpolation commutes with the exterior derivative (cf. [35,Section IV.27]). The exterior derivative of a k-Whitney form µ h ∈ W k h can be rewritten in terms of the incidence matrix d k ∈ R N k+1 ×N k (the discrete representation of the exterior derivative [28,39] This implies that the incidence matrix explicitly appears in the inner product where D k 1 is defined in (65), and in the duality product where G k 1 is defined in (66).

Discrete primal and dual Stokes-Dirac structure
In this section the discrete representation of the primal Stokes-Dirac structure (52) and the dual Dirac structure (53) are presented. Each of the two will verify a power balance involving inner products. By combining primal and dual variables a discrete Stokes theorem involving only duality products is obtained. This means that proposed discretization mimics the power balance in Eq. (27) at the discrete level.

Discrete primal Stokes-Dirac structure
The discrete version of the primal Stokes-Dirac structure (52) is obtained by replacing the flows and efforts with their discrete counterparts: The boundary flow f q−1 ∂,h is obtained by considering the trace of the associated polynomial family f q−1 ∂,h ∈ tr V q−1 s,h . Proposition 3. The system (83) verifies the following power balance This proves Eq. (84).
Using the algebraic realization of the inner product (62), the exterior derivative (65) and the trace (68), the algebraic realization of the discrete primal Stokes-Dirac structure (83) is expressed as where B q−1 s corresponds to a boundary control matrix, defined by

Discrete dual Stokes-Dirac structure
The discrete counterpart of the dual Stokes-Dirac structure (53) is given by In a reciprocal way with respect to the primal system, the boundary flow e p−1 ∂,h belongs to the trace of the p − 1 polynomial space, i.e. e p−1 ∂,h ∈ tr V p−1 s,h .
Proposition 4. The system (90) verifies the following power balance identity Summing Eqs. (92) and (93), it is found The algebraic realization of the discrete Stokes-Dirac structure (90) is expressed as The Whitney forms are of particular interest as they unveil the metrical and topological structure of the equations. Notice that the first equation of (97) and the second equation of (98) are purely topological, as they mass matrix can be factored out.

5.2.3.
Combining the primal and dual system: discrete power balance based on duality products The first equation in (83) and the second equation in (90) are a mere projection of the canonical Stokes-Dirac structure (18), i.e. the differential operator is enforced strongly. This means that by combining the primal and dual system a mimetic representation of the power balance in Eq. (27), which is based on duality products. This relies on the fact that finite element differential forms constitute a subcomplex of the de Rham complex. Proposition 5. The following discrete power balance, that combines variables from the primal and dual Stokes-Dirac structures holds Proof. Consider the first equation in (83) and the second equation in (90) The trimmed polynomial family forms a subcomplex of the de Rham complex [21,Lemma 3.8]. This assures that d e p−1 2,h ∈ V p s,h and de q−1 1,h ∈ V q s,h . Since the inner product is bilinear, it holds Since the inner product is non-degenerate, one has Taking the duality product with e q−1 1,h and e p−1 2,h and the two equations are summed, one obtains Given the fact that this holds ∀e q−1
This identity corresponds to the algebraic Stokes theorem reported in Eq. (75). In the Whitney forms case, the algebraic Stokes theorem in Eq. (75) reads Remark 7 (Connection with the construction of [26]). The relation (103) corresponds to the one reported in Proposition 4 in [26]. By factorizing the L matrices as , one can obtain an image representation of the Dirac structure in the minimal bond space (Proposition 5 in [26]). From the image representation an explicit port-Hamiltonian system can be constructed. The Hodge operator is subsequently built using the projection matrix obtained from the same factorization, which allows representing the constitutive equation. Here we stick to a FEEC construction, very much in the same spirit as [31], and by the simultaneous employment of primal and adjoint system the need of using a discrete Hodge is avoided altogether.

Discrete port-Hamiltonian systems
In this section the discrete representation of the primal and dual port-Hamiltonian systems are detailed.
Given a matrix A ∈ R n R ×n C with n R rows and n C columns, the notation [A] C R (with R = {r 1 , . . . , r #R }, C = {c 1 , . . . , c #C } index sets) indicates the matrix containing the rows and the columns of matrix A associated with the index set R and C respectively. Given a vector x ∈ R n R the notation [x] R indicates the vector containing the rows of vector x associated with the index set R only. In particular in the following sections, let us denote with a slight abuse of notation the index sets corresponding to the degrees of freedom of the interior of te domain and the subpartitions Γ 1 and Γ 2 respectively by I, Γ 1 , Γ 2 (cf. Fig. 2). With this notation, the degrees of freedom for the inputs defined in Eq. (25) are expressed as follows The degrees of freedom of the outputs, defined in (26), are given by . For sake of notational lightness, we stick to this slightly abusive notation. Concerning the actual implementation, it is crucial that the index set associated for each of the boundary variables satisfy Γ 1 ∩ Γ 2 = ∅. This requires that the geometrical entities at the intersection of the two are attributed to one or the other in an exclusive  Example of a 2D domain with splitting of the boundary due to the boundary conditions. In the case p = 2, q = 1 the duality is between a 0-form and a 1-form. Notice that the points at the intersection of Γ 1 and Γ 2 are attributed to Γ 1 as those will correspond to the essential boundary conditions for the 0-form.

Discrete primal port-Hamiltonian system
Starting from the weak formulation (56), the discrete primal port-Hamiltonian system is ob- The weak formulation (107) leads to the following discrete system

This system possesses an associated discrete Hamiltonian
Notice that the matrix [Ψ q−1 s ] Γ1 Γ1 is rectangular. Given definitions (112), (106), the statement is proven.
Remark 9 (The collocated nature of outputs y p−1 in a pHDAE setting). The collocated output in Eqs. (112) is a consequence of the strong imposition of the boundary conditions by direct assignment and represents the reaction to be applied at the boundary to guarantee to that the state follows the input trajectory. To highlight their collocated nature it is necessary to consider a differentialalgebraic framework and in particular port-Hamiltonian descriptor systems [46]. System (109) can be equivalently rewritten in the following canonical pHDAE form [46] Eẋ = Jx + Bu, where the state, input and output variables read The state now includes a Lagrange multiplier λ p−1 to enforce the constraint. The E, J, B matrices are given by (115) The variation of the energy is given by Notice that the output y q−1 2 corresponds to

Discrete dual port-Hamiltonian system
Starting from the weak formulation (58), the discrete dual port-Hamiltonian system is obtained: where u p−1 2,h ∈ tr V p−1 s,h | Γ2 and V q−1 s,h (Γ 1 ) is a polynomial space with boundary conditions The algebraic realization of the discrete weak formulation (117) reads The associated discrete Hamiltonian for this system reads (120) Proposition 7. The energy rate of the Hamiltonian in Eq. (120) (arising from the mixed discretization) is given by where y q−1 ] Γ2 is defined in Eq. and (105) Proof. From the dynamics (119), it is obtained Given definitions (122), (105), the statement is proven.
As in Remark 9, the strong imposition in System (117) of the boundary condition leads to a descriptor system with power balance where y q corresponds to

Recovering the power balance
For the case of uniform boundary conditions, it is sufficient to consider only one of the two mixed discretization Eqs. (109) and (119) to obtain the preservation of the discrete power balance (as in the Partitioned Finite Element method [31]). For example in the case in which Γ 1 = ∂M, Γ 2 = ∅, then the energy rate reads Conversely, if Γ 1 = ∅, Γ 2 = ∂M , then the energy rate reads However, in this general case of mixed boundary conditions, both systems (109) and (119) are needed in order to recover a discrete Stokes-Dirac structure and its associated power balance. The right hand side of the continuous power balance (21) written in terms of efforts can be expressed as P = e q−1 Its discrete version is then given by The dual field formulation is such to preserve the power balance at the discrete level, under some regularity assumption for the physical coefficients. This is due to the fact that each mixed system (107) and (117) contains an equation (the one belonging to the Stokes-Dirac structure) that is strongly verified, leading to a discrete conservation law.

Proposition 8.
If the tensors C p , E q verify the following regularity assumption then the discrete power satisfies Proof. Since the tensors satisfy the regularity assumption, the subcomplex property [21, Lemma 3.8] to the first line of (107) and the second line of (117) provides Taking the duality product against e q−1 1,h and e p−1 2,h provides The application of the discrete Stokes Theorem (74) then gives leading to the proof of the statement.

Time discretization
To ensure discrete conservation of energy, implicit Runge-Kutta methods based on Gauss-Legendre collocation points can be used [49]. These methods are also the only collocation schemes that lead to an exact discrete energy balance in the linear case [50]. The implicit midpoint method is here used to illustrate the time discretization.
Consider a total simulation time T end and a equidistant splitting given by the time step ∆t = T end /N t , where N t is the total number of simulation instants. The evaluation of a generic variable x at the time instant t n = n∆t is denoted by x n . The midpoint method is the simplest method in the class of collocation methods known as Gauss-Legendre methods. It applies to systems of the formẋ = f (x, t), This method is a one stage method, i.e. only uses the information at the previous time, and seeks for the solution of the implicit equation where x n+1/2 := xn+1+xn 2 and t n+1/2 := t n + ∆t 2 . The application of the implicit midpoint method to system (109) with direct assignment of the boundary degrees of freedom leads to the algebraic system The matrix A p bc corresponds to The b p bc vector incorporates the forcing due to the previous time step and the boundary conditions The implicit midpoint method, once applied to (119), leads to the system with and Proposition 9. Given the recursions implemented by systems (132) and (135), the following time discrete energy rate holds where y q−1 n+1/2 and y p−1 n+1/2 correspond a midpoint discretization of Eqs.
Proof. The time discrete energy rates are obtained by vector multiplication of Systems (132)  respectively. This computation provides the implicit midpoint discretization of Propositions 6 and 7.
The algebraic systems (132) and (135) satisfy a time discrete power balance, together with time discrete energy rates.
Proposition 10. Given the recursions implemented by systems (132) and (135), the regularity assumption reported in Eq. (128), and the following definition of time discrete power Proof. The time discrete power balance corresponds to the implicit midpoint time discretization of Proposition 8.

Numerical experiments
In this section, the dual field discretization methodology proposed in this work is tested for the wave and Maxwell equations in a box-shaped three-dimensional domain The boundary sub-partitions are selected to be Since we want to rely on existing and well-established librairies like FEniCS and Firedrake, the equations need to be translate into vector calculus operations. By introducing the musical isomorphism, given by the flat and the sharp operator (cf. Appendix B), and the isomorphism β converting vector fields in n − 1 forms, the commuting diagram in Fig. 3, that provides the link between the de Rham complex and the standard operators and Sobolev space from vector calculus, is obtained. Since the manifold is a subset of the Euclidean space, the metric tensor is the identity g ij = δ ij . The finite element arising from the trimmed polynomial family on the computational mesh T h are then equivalent to the well known continuous Galerkin (or Lagrange) elements P − , as illustrated in Figure 4.
s,h grad curl div Figure 4: Equivalence between finite element differential forms and classical elements.
The finite element library Firedrake [37] is used for the numerical investigation.
Remark 10. As argued in [45], the method of manufactured solution is not suited for boundary control problems. This is due to the fact that for finer discretizations the volume terms due to the forcing will dominate over the boundary terms. This is the reason why the numerical tests are set up considering an eigensolution rather than a manufactured one induced by a forcing.

The acoustic wave equation in 3D
The acoustic wave equation corresponds to the case p = 3 and q = 1. Using the same notation as in Sec. 3.1, the energy variables are the top-form v 3 := α 3 (corresponding to the pressure) and the one-form σ 1 := β 1 (corresponding to the linear momentum). If the physical coefficients are normalized to one (this can be easily achieved by re-scaling the time with respect to the speed of sound), the Hamiltonian is given by with its variational derivatives given by leading to the pH system Given the functions g(x, y, z) = cos(x) sin(y) sin(z), an exact solution of (143) is given by The exact solution provides the appropriate inputs to be fed into the system The employment of the dual field discretization leads to the resolution of two systems: • the primal system (109) of outer oriented variables v 3 and σ 2 ; • the dual system (119) of inner oriented variables v 0 and σ 1 ; Each variables is discretized using the associated finite element differential forms: • discontinuous Galerkin elements DG s−1 for v 3 h ; • Raviart Thomas elements RT s for σ 2 h . • continuous Galerkin elements CG s for v 0 h ; • Nédélec elements of the first kind NED 1 s for σ 1 h ;

Energy Conservation properties
First the conservation properties of the scheme are verified against the exact solution (145). The test is performed using N el = 4 elements for each side of the box-shaped domain and polynomial degree s = 3. Concerning the time discretization, the total simulation time T end = 5 and the time step is taken to be ∆t = T end 200 .
The energy rate conservation, expected by Prop. 9, is reported in Fig. 5. Fig. 6a shows the fulfilment of the discrete power balance (Prop. 10), whereas Fig. 6b provides the error between the discrete and exact power flow. The power flow error due to the polynomial interpolation is lower than 10 −4 . For what concerns the energy behaviour, three different energies are considered The first two energies are associated with the mixed discretization (109) and (119) respectively, whereas the last one is constructed using the dual field formulation. The errors on the energy rate and variation of energy, shown in Fig. 7, are of the order 10 −5 , assessing the performance of the dual field formulation. In particular, one can notice that the energy H T 2 stays between H 3 h , H 1 h at all times. The variation of energy is also computed using the power flow As shown in Fig. 7b this quantity is less affected by the error for this particular test case.

Convergence results
In this section, the convergence rate of the different variables is verified against the analytical solution (145). The error is evaluated in the L 2 norm at the ending time. Concerning the time discretization, the total simulation time T end = 1 and the time step is taken to be ∆t = T end 100 .
In Fig. 8 the L 2 error trend against the exact solution is reported for all variables. It can be noticed that the error goes as h s , exception made for v 0 (cf. Fig. 8b) that for s = 1, 2 stays between h s and h s+1 . Indeed, the presence of mixed and time-varying boundary conditions leads to a lower convergence rate than the expected theoretical order of h s+1 for homogeneous boundary conditions (see [51] for the error analysis of the Hodge wave equation under one case of homogeneous conditions). The other variables exhibit the same convergence trend as predicted by the analysis in [51]. The L 2 difference between the dual representation of the variables is reported in Fig. 9. The difference of the dual representation converges as h s . This is in accordance with the results obtained in [33], where the dual field formulation is employed to solve the Navier-Stokes equations in periodic domains.

The Maxwell equations in 3D
The Maxwell equations corresponds to the case p = 2, q = 2. The energy variables correspond to the electric displacement two form d 2 = α 2 and the magnetic field b 2 = β 2 . The Hamiltonian reads where ε ∈ R is the electric permittivity and µ ∈ R is the magnetic permeability. The variational derivative of the Hamiltonian are given by Variables e 1 , h 1 are the electric field and the magnetizing field respectively. Since the reduction of the constitutive equation is such to keep only the efforts variables and their duals, the following dynamical system is obtained.      where e 2 = e 1 , h 2 = h 1 . Given the functions where ω = √ 3c and c = ( √ µε) −1 , the system (149) is solved by the eigenmode The exact solution provides the appropriate inputs to be fed into the system The employment of the dual field discretization leads to the resolution of two systems: • the primal system (109)

Conservation properties
The conservation properties of the scheme are verified against the exact solution (151). The test is performed using N el = 4 elements for each side of the box-shaped domain and polynomial degree s = 3. Once again, the total simulation time T end = 5 and the time step is taken to be ∆t = T end 200 .
An important feature of the Maxwell equations (149) is that they verify the following constraints d 2 e 2 = 0, This result follows by taking the exterior derivative of each line of system (149), under the assumption that the initial conditions respect these constraints. Mixed finite element strategies, like the ones proposed in [44,52,53], cannot satisfy both constraints as they do not employ a dual representation for each variable. Instead, the dual field formulation naturally capture this aspect as shown in Fig. 10. The energy rate conservation and discrete power balance are reported in Figs. 11 and 12a respectively. The numerical test confirms once again the expected behaviour. The power flow error due to the polynomial interpolation is of the order of 10 −5 (cf. Fig. 6b).
For what concerns the energy behaviour, three different energies are once again considered Fig. 13 shows that the error on the energies are of the order 10 −5 . The dual field energy H T 2 stays in the middle of H 2 h , H 2 h . The variation of energy is also computed using the power flow Indeed this variation of the energy is not the most precise (cf. Fig. 13b). A rigorous error analysis is needed to assess the conditions under which one of these energies perform better. All variables converge with a rate given by h s (see Fig. 14). However, it appears that for s = 2, 3 variables e 2 h and h 2 h (Figs. 14a, 14c) the convergence rate is a little less than h s . A rigorous error analysis is needed to gain more insight about the observed behaviour. An a priori analysis using a mixed finite element scheme can be found in [44]. However, therein only homogeneous and uniform boundary conditions are considered.

Additional insights on the choice of the dual variables
Finally, we conclude by a short discussion on an alternative way of defining the diffeomorphism Φ in (30) which relates the original energy variables to the dual ones and gives rise to the adjoint Stokes-Dirac structure. Differently than presented in Sec. 3.4, we herewith take a definition of the Adjoint Dirac-Structure which will create an interesting symmetry between the primary system and the adjoint system. This will have as a consequence that the material operators A p , B q will then appear in the adjoint Dirac-Structure. This is not useful for discretization purposes, but is much more natural from a physical point of view, because the metrical property of space are strictly related not only to the Hodge, but to the coupling of the Hodge with the material properties.
We herewith then redefine the dual states (α n−p , β n−q ) to achieve a redundant representation of the state space in such a way that the dependencies on both the metric and material properties disappear. We will also redefine the adjoint Hamiltonian H(α n−p , β n−q ) and we will do it here in such a way that the sum of the original and adjoint Hamiltonian is independent from the metric and material properties.    Proposition 11. Given any k-form ω k and a symmetric positive definite isomorphism A k between k forms, if we define ω n−k := A k ω k , and C n−k := (−1) k(n−k) A −1 (which is symmetric positive definite), we obtain the following relations: and (A k ω k , ω k ) M + ( ω n−k , C n−k ω n−k ) M = 2 ω k | ω n−k M .
Proof. Considering that A k is an isomorphism, we have that ω k = (−1) k(n−k) A −1 ω n−k . Since the Hodge star is an isometry, it is found Furthermore, one has ( Aω k , ω k ) M = ω k | ω n−k M , leading to the second equation.
We can therefore use the previous results by defining the following flow relations: where A n−p = (−1) p(n−p) A p , B n−q = (−1) q(n−q) A q . If A p , A q are not regular enough, then they cannot be differentiated and it is necessary to invert them and bring them to the left side (158) This system can be then put into weak form considering the integration by parts applied to the codifferential. Relation (153) directly gives the representation which can be used to define the Hamiltonian of the adjoint system, which then represents what in physical system theory is called the co-energy H(α n−p , β n−q ) = M 1 2 C n−p α n−p ∧ α n−p + 1 2 B n−q β n−q ∧ β n−q .
Furthermore (154) immediately shows that by creating a double representation of the system with a different dual state, the sum of the energy and co-energy, can be expressed as purely as function of the states and no extra metrical properties H( α p , β q ) + H(α n−p , β n−q ) = M α p ∧ α n−p + β q ∧ β n−q , which shows that the sum of the energy and co-energy is independent of the metric properties expressed by A p and A q achieving a perfect symmetry. It can also be seen that in this alternative definition of the adjoint Dirac-Structure, the metric properties of space represented by the Hodge are always taken together with the physical properties of space represented by the A operators as it would be expected from a physical point of view. The insights presented above could be instructive in extending our proposed discretization scheme to nonlinear port-Hamiltonian system.

Conclusion
In this contribution, the dual field formulation is employed for the systematic discretization of linear port-Hamiltonian systems under generic boundary conditions. The proposed methodology is entirely based on the finite element exterior calculus framework. The dual field formulation solves the problems associated with the construction of a discrete Hodge operator (that typically requires dual topological meshes to preserve its isomorphic character) by relying on the adjoint system. The employment of the adjoint system introduces the boundary conditions explicitly by means of the integration by parts formula. This leads to two decoupled mixed discretizations that, once solved, allow retrieving a discrete power balance, regardless of the underlying boundary conditions. This guarantees that the proposed discretization method gives rise to a Dirac structure. This is of crucial importance for multiphysics applications, as well as the design of model-based control laws.
This methodology opens the door to a number of interesting developments. As argued in [33], the employment of dual representations for the unknowns provides useful indicators for adaptive meshing, as one can use the norm of the difference between dual variables as an a posteriori estimator. Furthermore, by using a time staggered discretization as in [33] the boundary conditions could be imposed in a completely weak manner. This would allow the construction of an explicit state-space model, thus avoiding the complications associated with differential algebraic systems. An interesting aspect concerns the usage of algebraic dual polynomials (proposed in [54]) so that dual solutions can be represented in a pair of algebraically dual polynomial spaces.
The dual field formulation has been successfully employed to tackle the rotational term of the Navier-Stokes equations in a linear manner, leading to a computationally efficient scheme that conserves mass, energy and helicity. For this reason, we expect a non linear extension of this method to be feasible and competitive with respect to state of the art solutions. An interesting development concerns the extension of the proposed methodology to elasticity problems. These problems require a non trivial extension of the canonical Stokes-Dirac structure, that is based on the de Rham complex, as the differential operators included in the underlying complex, the elasticity complex, are not topological but metrical.