Theory and Implementation of Coupled Port-Hamiltonian Continuum and Lumped Parameter Models

A continuous Galerkin finite element method that allows mixed boundary conditions without the need for Lagrange multipliers or user-defined parameters is developed. A mixed coupling of Lagrange and Raviart-Thomas basis functions are used. The method is proven to have a Hamiltonian-conserving spatial discretisation and a symplectic time discretisation. The energy residual is therefore guaranteed to be bounded for general problems and exactly conserved for linear problems. The linear 2D wave equation is discretised and modelled by making use of a port-Hamiltonian framework. This model is verified against an analytic solution and shown to have standard order of convergence for the temporal and spatial discretisation. The error growth over time is shown to grow linearly for this symplectic method, which agrees with theoretical results. A modal analysis is performed which verifies that the eigenvalues of the model accurately converge to the exact eigenvalues, as the mesh is refined. The port-Hamiltonian framework allows boundary coupling with bond-graph or, more generally, lumped parameter models, therefore unifying the two fields of lumped parameter modelling and continuum modelling of Hamiltonian systems. The wave domain discretisation is shown to be equivalent to a coupling of canonical port-Hamiltonian forms. This feature allows the model to have mixed boundary conditions as well as to have mixed causality interconnections with other port-Hamiltonian models. A model of the 2D wave equation is coupled, in a monolithic manner, with a lumped parameter model of an electromechanical linear actuator. The combined model is also verified to conserve energy exactly.


Introduction
As computational power increases and the desire for ever more complex models grows, the need to have an energy-conserving mathematical framework for multiphysics, multidomain problems is becoming apparent. Many fields require complex couplings between continuum and lumped parameter models (LPMs) including physiology, aerospace, vehicle dynamics, robotics, and many more. The coupling of these models is typically done in an iterative manner, which has multiple disadvantages. Two major downsides are the failure of iterative couplings to conserve energy [35] and the difficulty in finding a stable combination of timesteps for the coupled models [45]. Another equally large disadvantage of the iterative approach is the difficulty of implementing control algorithms without a monolithic statespace matrix. A fully coupled, monolithic approach that conserves the structural properties of the lumped parameter and continuum models is desirable.
In this paper we discuss the energy conservation of a continuum system described by a partial differential equation (PDE) coupled with a LPM. Specifically, a model is derived for the linear 2D wave equation with an electromechanical (EM) linear actuator driving one of the boundaries. This example also shows that the port-Hamiltonian (PH) formalism can be used to interconnect various types of models while conserving energy flow through ports/boundaries. PH is a proven method for the modelling and control of complex multiphysics systems. Over the last 20 years, a significant amount of work has been done on infinitedimensional PH methods for the modelling of continuum systems [17,42]. Some PDEs that have been modelled using the PH framework are transmission line equations, shallow water equations and Timoshenko beam equations [14]. More recently, both 2D and 3D models have been implemented using the PH framework [38,41,44].
In Sect. 2 the spatial discretisation is proven to conserve the Hamiltonian in the same manner as the exact equations. Section 3 formulates the discretisation of the temporal domain with well known symplectic integrators: the symplectic Euler (SE) method, the symplectic (implicit) midpoint (SM) method and the Störmer-Verlet (SV) (leapfrog) method. A symplectic method conserves volume in phase space, which results in bounded conservation of the Hamiltonian [20]. The combined spatial-temporal discretisation, therefore, conserves the Hamiltonian structure of the governing equations. A similar class of methods that conserves the Hamiltonian structure is the class of multi-symplectic methods. Multiple groups have applied multi-symplectic methods to both the linear and non-linear wave equations. Reich used Runge-Kutta finite difference schemes in both space and time [34], McLachlin compared multiple methods, including spectral methods [28], and Brugnano used a Hamilton Boundary Value Method [5]. McDonald shows that multi-symplecticity preserves the travelling waves of hyperbolic equations [27]. The Partitioned Finite Element Method (PFEM) in the paper by Cardoso-Ribeiro et al. [11] combines a finite element spatial discretisation with an SV time integration scheme in a way that conserves energy but requires Lagrange multipliers to implement mixed boundary conditions [21]. Brugnoli et al. successfully apply this approach to Mindlin and Kirchhoff plate models [7,8]. Brugnoli et al. [9] have also introduced another method for applying mixed boundary conditions to PFEM. Their method requires discretisation of the spatial domain into separate sections, each section having one type of boundary condition. The Stokes-Dirac structure [43] and interconnection property of PH methods is then used to combine the multiple sections, creating a full system with mixed boundary conditions. Similar to our method, Kotyczka uses a finite element discretisation and symplectic time integration in a way that conserves the Hamiltonian structure and allows for mixed boundary conditions [22], however, user-defined parameters in the method must be tuned for accurate results. The method introduced in this paper combines desirable attributes of Cardoso-Ribeiro's, Brugnoli's, and Kotyczka's methods and provides a Hamiltonian-conserving, symplectic method that allows for easily implemented mixed boundary conditions, port-based boundary coupling, and does not require tuning of userdefined parameters.
To ensure conservation of energy flow through the boundaries, a weak boundary condition implementation is used for the Dirichlet conditions, similar to the way Neumann conditions are typically implemented in finite element methods. Weak boundary conditions are implemented in the variational form and provide many benefits for a finite element formulation. One of these benefits is that it simplifies the implementation by not having to directly prescribe the degrees of freedom (DOFs) at the boundary. This can have benefits in applying multiple types of boundary conditions, including no-slip conditions for the Navier-Stokes equations [3]. Also, no manipulation of the solution matrix is required to prescribe the DOFs. The ability to have mixed boundary conditions is also extended to allow mixed causality interconnections with other PH models. By showing that our spatial discretisation is equivalent to a coupling of canonical PH models with either Neumann or Dirichlet boundaries and by following the interconnection methods of Brugnoli [6], we calculate the causal boundary connections with a LPM. The ability to use the canonical forms to calculate the power conserving interconnection allows mixed causality boundary interconnection between other PH models.
Continuous Galerkin and hybridizable discontinuous Galerkin methods couple well with PH methods due to the ease of ensuring structural properties of the models and the simplicity of coupling different models at the boundary. Brezzi and Fortin give a good overview of Galerkin methods [4]. Cockburn developed a unifying framework for Galerkin methods [12] that multiple authors have extended. McLachlan extends Cockburn's discontinuous Galerkin (DG) methods into a formulation that allows proof of multi-symplecticity for elliptic equations [29]. Sanchez uses a hybridizable discontinuous Galerkin approach to obtain a wave equation discretisation similar to our method, that is proven to be Hamiltonian-conserving and symplectic in time [36]. We prove the same qualities for our method, a continuous Galerkin approach that has the important added novelty of being able to be coupled with arbitrary PH LPMs. The conservativity and modular approach of this method is thus ideal for a wide range of real-world problems.
Energy conservation in the Galerkin method is extremely important for real-world applications that require long-time simulations, such as in geophysical fluid dynamics. Bauer uses a Poisson bracket approach to prove conservation of energy for a Galerkin discretisation of the rotating shallow water equations [2]. Extending on Bauer's work, Eldred uses the Galerkin method coupled with a Poisson time integrator to conserve energy when modelling the thermal shallow water equations [15]. Both of these models take advantage of Hamiltonian-conservation to give good long-time prediction for geophysical fluid dynamics applications.
Modelling of wave propagation also has applications in multiple biological fields. In this paper, the wave equation models homogeneous, linear-elastic media. For extensions with heterogeneous materials, typical of most biological materials, see the paper by Serhani [38]. Elastography [31] is an interesting application that requires a type of inverse modelling to identify the elasticity of heart tissue. A coupled LPM-continuum model could improve the current elastography methods. There has also been work done on non-linear models for wave propagation through biological tissue [13].
The model in this paper is implemented with the software FEniCS [1], which is a tool for automated scientific computing that focuses on solving PDEs. FEniCS allows for a very efficient implementation of finite element methods specified in a weak form. Another strong attribute of FEniCS is its automatic differentiation, which is valuable for inverse problems and for control.
The section outline of the paper is as follows. Section 2 details the PH form of the wave equation, proves that the discretisation conserves energy, and proves that the discretisation retains the Stokes-Dirac structure. The FEniCS implementation and the results showing energy conservation are shown in Sect. 3. Section 4 validates the wave equation against an analytical solution by showing spatial, temporal, and eigenvalue convergence. An EM model is introduced in Sect. 5. In Sect. 6 the wave equation is coupled with the EM model in a monolithic approach that conserves the canonical PH structure. Finally, in Sect. 7, the results of the combined model are shown and the energy conservation of the model is discussed. Appendices A and C detail the Python code for Sects. 3 and 6, respectively.

The Wave Equation
This section details the weak form PH discretisation of the wave equation with constant propagation speed, the conservation of energy proof, and the proof that the discretisation ensures a Stokes-Dirac structure. The 2D wave equation in Cartesian coordinates is used throughout this paper, however, all results in this section naturally extend to the 3D wave equation. The basic linear wave equation is where w is the wave amplitude, x denotes the spatial coordinates, and t is the time. To model the simplified wave propagation in an elastic membrane we define c 2 = k w /ρ w as the wave speed squared, which is a constant function of the material density (ρ w ) and stiffness (k w ). To transform Equation (1) into PH form, state variables,p, the momentum, andq, the strain, are chosen asp This transforms the second-order equation into a system of (n + 1) first-order equations, where n is the number of spatial dimensions. Note that the tilde overscript is used to denote exact variables, to distinguish them from the approximate functions used in subsequent sections. Using the PH notation of flow and effort variables, the time derivative of the state variables,f p andf q , are defined as flows The Hamiltonian functional and the Hamiltonian density are respectively given by where is an open, bounded spatial domain with a Lipschitz-continuous boundary, ∂ . The effort variablesẽ p , the velocity, andẽ q , the stress, are defined as the variational derivatives of the Hamiltonian density, For a functional that only depends on its states, not on their spatial derivatives, the variational derivatives are equal to the partial derivatives of the integrand. Transforming Equations (1) and (3) to (6) into a PH structure gives where the div and grad operators make up the formally skew-adjoint J operator. For a proof of the skew-adjointness of J see the work by Trenchant et al. [40].
Theorem 2.0.1 Equation (7) is energy-conserving. i.e., the rate of change of the Hamiltonian is equal to the energy flow through the domain boundary.
Equation (8) satisfies the well known bond-graph and PH condition, that the product of the effort and flow variables equals the power [42]. Substituting Equation (7) into Equation (8) and using integration by parts giveṡ Therefore, the rate of change of the Hamiltonian is only dependent on the boundary terms and thus the Hamiltonian is conserved within the internal domain.
respectively, or for the opposite causality the inputs and outputs are respectively,

Weak Form
In this section, the discretised weak form of the wave equation is derived. First, we define the L 2 inner product over the domain, , and the boundary, ∂ , as Using the same function spaces as in Cardoso-Ribeiro's work [11], approximate flow and effort functions are introduced, f p , e p ∈ H 1 ( ) , f q , e q ∈ H div ( ) . (12) In the succeeding sections, we show that this choice of function spaces allows us to combine important features of Cardoso-Ribeiro's and Kotyczka's [22] methods. Cardoso-Ribeiro's methods ease of implementation is combined with Kotyczka's methods ability to implement mixed boundary conditions without Lagrange multipliers. Substituting these approximate functions for the exact flows and efforts in Equation (7) and taking the inner product with the test functions The right-hand side of Equations (13a) and (13b) are then integrated by parts to give These equations are now in a form where the Galerkin method can be applied. To do this, the approximate flow and effort functions in Equation (12) are defined from the discrete flow (f p ∈ R Np ,f q ∈ R Nq ) and effort (ê p ∈ R Np ,ê q ∈ R Nq ) vectors and the vectors of globally defined basis functions (ϕ p , ϕ q ) as shown here, where N p and N q are the number of DOFs stored by the discrete vectors. A hat over a variable is used to denote the vector of discrete values, i.e.,ê p is the column vector of discrete DOF values for the scalar field e p . To make the method Galerkin, v p and v q are discretised with the same basis functions as f p and f q , respectively. The basis function families that we use are Lagrange for ϕ p and Raviart-Thomas [33] for ϕ q , however, any basis functions that satisfy the function spaces in Equation (12) are suitable. In Equations (15) and (16), the sizes of ϕ p and ϕ q are N p × 1 and N q × 2, respectively. The notation, ·, · is used for the standard inner product on R, as opposed to ·|· , the L 2 inner product. For the lowest order Lagrange and Raviart-Thomas elements,ê p andf p are stored at nodes andê q andf q are stored at edges. For details on higher order elements see the FEniCS book [26]. Substituting the approximate functions from Equations (15) and (16) as well as the corresponding test functions into Equation (14a) gives where ψ p = ϕ p | ∂ represents ϕ p evaluated at the boundary and ψ q = ϕ q · n| ∂ represents ϕ q evaluated in the normal direction at the boundary. More clearly, the basis functions for the boundary terms satisfy the following relations, Substituting the approximate functions from Equations (15), (16) and (18), as well as the corresponding test functions into Equation (14b) gives The matrices in Equations (17) and (19) are given as where the inner product acts elementwise for the matrix as Applying a typical weak form approach, Equations (17) and (19) must hold for anyv p ∈ R Np andv q ∈ R Nq . Therefore, according to the fundamental theorem of variational calculus, the following matrix system of equations holds, M pf p = K pêq + L pêq , When implementing Equation (22), it is essential to include the discrete dynamic and constitutive laws which coincide with Equations (3) and (6), respectively. The dynamic and constitutive laws are respectively, where I is the identity matrix. Also,p andq are the discrete vectors of momentum DOFs and strain DOFs, respectively. This formulation assumes constant material properties and therefore a constant Q matrix. Finally, substituting Equations (23) and (24) into Equation (22) we get the matrix system of equations, In Sect. 2.3 we prove that Equation (25) can be written in two different (one for Dirichlet and one for Neumann boundary conditions) canonical PH forms and that it represents a non-degenerate Stokes-Dirac structure.

Discrete Conservation of Power Proof
To proceed with a conservation of power proof we first define the approximate Hamiltonian, H , and Hamiltonian density,Ĥ in the same way as the exact functions that were defined in Equations (4) and (5), The variables p and q have the same basis functions as e p and e q respectively. Substituting in the discrete vectors,p andq and their corresponding basis functions giveŝ The relationship between the approximate state variable functions and the approximate effort functions from Equation (16) can be found by taking the partial derivative ofĤ(p, q) with respect to p and q.
The corresponding relation between discrete effort variables and state vectors iŝ To retain the structure of the continuous system, the discretised equations must have the same energy-conserving structure as the continuous equations. Therefore, the rate of change of the Hamiltonian must only depend on the boundary variables, as in Theorem 2.0.1. The following conservation of energy proof is influenced by the thesis of Kotyczka [22] and the paper by Cardoso-Ribeiro [11]. To prove Theorem 2.2.1 a mapping between general variables and boundary variables must be formulated. Following the work of [22], we decompose L p into T p , S q and L q into T q , S p , as shown, The matrix T p is simply a mapping from allê p DOFs of the mesh to the DOFs that havê e p defined at the boundary. Similarly, T q is a mapping from allê q DOFs of the mesh to the DOFs that haveê q defined at the boundary. Both T p and T q are made of zeroes and ones and are semi-orthogonal, therefore T p T T p = I Np and T q T q T = I Nq . For first order Lagrange and Raviart-Thomas elements, the identity matrix I Np has a number of rows and columns equal to the number of nodes (N p ) and I Nq has a number of rows and columns equal to the number of edges (N q ). Calculation of S p and S q is done trivially by using the semi-orthogonality of T p/q and Equation (30), giving S q = T p L p and S p = T q L q . The matrix S p is a mapping from momentum effortsê p to boundary strain flowsf qb . Similarly, S q is a mapping from stress effortsê q to boundary velocity flowsf pb . The matrix mappings of T p/q and S p/q are summarised asê pb = T pêp ,ê qb = T qêq , To prove Theorem 2.2.1 we need a small lemma.

Lemma 2.2.2 (Lemma 1 [24])
Applying integration by parts in reverse then gives therefore, Proof of Theorem 2.2. 1 We start by taking the time differential of the Hamiltonian in the usual way.Ḣ The approximate efforts and flows from Equations (15) and (16) are then substituted for the effort and flow functions to give the discrete Hamiltonian rate of changė Equation (22) is then used to givė Lemma 2.2.2 can now be used. Following that, the definitions of the boundary flows and efforts in Equation (31) can be used to complete the first half of the proof, and identically for the second half of the proof,

Corollary 2.2.3
The rate of change of the Hamiltonian of the discrete system conserves energy in the same way as in Theorem 2.0.1, for the continuous system.

Proof
Beginning with Equation (33) and substituting in L p from Equation (20) giveṡ Substituting Equation (18) giveṡ which confirms that the rate of change of the discrete system Hamiltonian is a function of the efforts at the boundary. This correctly coincides with Theorem 2.0.1.
The approximate system inputs and outputs are thus defined in the same way as Equations (9) and (10), either or

Stokes-Dirac Structure and the Canonical Port-Hamiltonian Form
In Equation (22) we introduce the model in a non-canonical PH form, however, in this section we prove that Equation (22) is equivalent to a structure-preserving coupling of canonical Port-Hamiltonian forms. This also proves that the discretisation ensures a Stokes-Dirac structure [43], therefore, conserving the structure of the continuous equations, Equation (7). To formulate the system in an input-state-output PH form, we first define the input and output functions at the boundary corresponding to Equations (34) and (35) in terms of their discrete vectors (û q ,ŷ p ,û p ,ŷ q ) and basis functions (θ q , θ p ).
where θ p and θ q contain the entries of ψ p and ψ q corresponding to the boundary DOFs. The first canonical PH form for Neumann boundary conditions can be set up by using which was formulated in Equation (32). Inserting the transpose of Equation (38) into Equation (22) gives We then multiply the output equation in Equation (34) by the v yq = v yq |θ q ∂ N trial function, where a subscript N denotes a Neumann boundary. Following this, we integrate over the boundary then use the fundamental theorem of variational calculus to get, where M yp = θ q |θ p ∂ N . Also, B qb is defined as, Combining Equations (39) to (41) gives the PH canonical form, ⎡ ⎢ ⎣ Including the dynamic and constitutive laws from Equations (23) and (24) gives the input- where M qb is the mass matrix, J qb is the skew-symmetric matrix, and Q qb is the constitutive law matrix. The infinite-dimensional canonical form corresponding to the discretised canonical form in Equation (42) is where, M qb is the mass operator and J qb is the skew-symmetric operator. Technically, Equations (43) and (44) are canonical forms for domains with Neumann conditions only. In the following we formulate the canonical form for Dirichlet conditions. We take the transpose of Equation (38) and use the fact that L q = L T p to give Equation (45) can be substituted into Equation (22) and the discretisation of y q can be done in the same way as for y p in Equation (40) to give the Dirichlet boundary condition equivalent of Equation (42). This discrete canonical form is ⎡ and the input-state-output PH form is where M yq = θ p |θ q ∂ D and a subscript D denotes a Dirichlet boundary. The canonical PH form in Equation (47) has a mass matrix, M pb , a skew-symmetric matrix, J pb , and a constitutive law matrix, Q pb . The interconnection matrix, B pb is defined as The infinite-dimensional canonical form corresponding to Equation (46) is where, M pb is the mass operator and J pb is the skew-symmetric operator.
Equations (47) and (49) are only canonical PH forms for boundaries with Dirichlet conditions. For a domain with mixed boundaries, our system in Equation (25) needs to be equivalent to a canonical PH system, or in our case a combination of canonical PH systems. Any closed domain with mixed Neumann and Dirichlet boundary conditions can be subdivided into subdomains with only Neumann or only Dirichlet boundary conditions. This idea has been taken advantage of in Brugnoli et al.'s work [9], where the authors numerically segment the domain and apply PFEM to each section. Here the idea is only used conceptually to prove that our system in Equation (25) is equivalent to a combination of canonical input-state-output PH formulations and therefore, by the compositionality property, retains the Stokes-Dirac structure of the analytic equations, Equation (7). This means that mixed boundary conditions can be implemented in the weak form, as detailed in Sect. 3, and that the resulting system is a Stokes-Dirac structure. It is also important to note that our formulation is non-degenerate, due to M p and M q being full rank. The matrices turn out to be full rank because we use the same basis functions for effort and flow functions. Proposition 1 of Kotyczka's paper [24], uses a similar compositionality argument to allow mixed boundary conditions and mixed causality at the boundaries. However, in Kotyczka's work, the nonfull-rank M p and M q matrices cause the requirement of user-defined parameters in-order to form a non-degenerate Stokes-Dirac structure. Proof First we subdivide the domain so that each connected Neumann boundary is in a domain denoted Ni that does not connect or overlap with any Nj (j = i) and is not connected to any other external boundaries. A simple example of a subdivision is shown in Fig. 1. Each of the i subdomains has inputs and outputs split up into u qi , y pi at the external boundary, and u int qi , y int pi at the internal boundary. The remainder of the domain is therefore connected and only has Dirichlet boundary conditions, we denote this subdomain D .
The inputs and outputs of D are split up into u pk , y qk at the external boundaries and u int pi , y int qi at the internal boundaries that border Ni . The canonical form of Equation (44) is used in the Ni domains and Equation (49) is used in the D domain. Note that the canonical forms are modified accordingly, to split the inputs and outputs into internal and external parts. The causal interconnection relations [6,24] at each internal boundary can then be written as This ensures conservation of energy flow between the subdomains due to the power conserving inner product, Therefore, since each subdomain has a PH Stokes-Dirac structure, the power conserving interconnection ensures the total system is also a Stokes-Dirac structure by compositionality. Also, each subdomain has the correct canonical form for its Neumann or Dirichlet boundary conditions. Lastly, any connected domain is equivalent to a decomposition in the way we have described. The combination of these results proves that Equation (25) is equivalent to a structure-preserving combination of canonical PH forms for any connected domain with mixed boundary conditions.
As will be seen in Sect. 6 the infinite-dimensional canonical form in Equation (49) can be used to determine the power conserving interconnection between other canonical PH systems, where the interconnection enforces a Dirichlet condition on the wave domain boundary. Similarly, the infinite-dimensional canonical form of Equation (44) can be used to find the power conserving interconnection for connections that assign a Neumann condition on the wave domain boundary. Due to the conclusion of Theorem 2.3.1, that any domain is equivalent to a subdivision of subdomains with either the canonical form of Equation (49) or Equation (44), both Dirichlet and Neumann interconnections can be implemented on a domain. This means that the system developed in this paper can have mixed boundary conditions as well as mixed causality interconnections between other PH systems.

Wave Equation FEniCS Implementation
In this section, the FEniCS implementation of Equations (14a) and (14b) on a rectangle domain is detailed. Appendix A supplements this section by detailing the Python code for the implementation. A schematic of the wave domain is shown in Fig. 2.
An unstructured, triangular mesh was created over the domain with FEniCS meshing software. The boundary conditions on the domain are set as where the inputs in Equations (34) and (35) for the left, right, and middle boundaries are defined respectively as Dirichlet conditions are applied to both the left and the right boundaries. The left boundary is set as an input condition whereas the right boundary is given a fixed zero value. The top and bottom boundaries have a zero-flux Neumann condition applied. Separating the boundary terms in Equations (14a) and (14b) for each boundary condition and reverting to integral rather than inner product notation gives In Equations (54a) and (54b), the terms inside the boundary integrals are evaluated at their respective boundary. Substituting the state variables f p = −ṗ, f q = −q, e p = p/ρ w , e q = k w q and the boundary terms from Equations (52a)-(52c) gives These equations can be implemented in FEniCS, which automatically generates a matrix system of equations in the form of Equation (25). However, to solve these equations the system must also be discretised in time. Symplectic time integration schemes conserve the symplectic structure of the continuous equations and approximately conserve the Hamiltonian, therefore, they are the natural choice for the temporal discretisation. The symplectic Euler (SE) time integration scheme is applied to Equations (55a) and (55b) to give where the m superscript denotes the variable at the previous time step. The L, R, and M subscripts denote variables at the left, right, and middle boundaries, respectively. The SE scheme combines an explicit step for Equation (56a) and an implicit step for Equation (56b). When the Hamiltonian is separable the SE scheme is semi-explicit, meaning Equation (56b) could be solved explicitly after the solution of Equation (56a). However, due to the ease of implementation in FEniCS, the equations are solved in one implicit step. By Theorem 2.2.1, Equations (55a) and (55b) conserve the Hamiltonian. Combining this with SE integration gives a discrete system that retains the Hamiltonian structure of the continuum equations and conserves energy for large times, as further discussed in Sect. 3.1.
The energy bound on symplectic methods is, in general, proportional to O(( t) r ), where r is the order of the time integration scheme [20]. To compare the energy bound between different order methods, a second-order symplectic scheme, the Störmer-Verlet (SV) method [19] is also implemented, Solving Equations (57a)-(57c) requires two system of equation solves per time step, Equation (57a) is solved to get q 1/2 then Equations (57b) and (57c) are solved for q and p. Although in the general case the energy residual is bounded for symplectic methods, since the example we are modelling is linear we can improve this result and get exact energy conservation. We do this with the symplectic midpoint (SM) method, which conserves all quadratic invariants for linear systems [20]. The variational form of the SM scheme is

Wave Results
In this section, the model that is detailed in Sect. 3 and proven to conserve energy by Theorem 2.2.1 has been implemented with a range of different time integration schemes. Our method, which has a weak Dirichlet boundary condition implementation is compared to an implementation with Dirichlet boundary conditions implemented in the typical strong form. This is done to show that the naive setting of boundary conditions in a strong manner is detrimental to energy conservation. The strong Dirichlet implementation modifies the matrix system of equations generated by FEniCS to directly enforce the boundary condition values. This differs from our weak boundary Dirichlet implementation, where we apply the boundary conditions by specifying the boundary integral terms in Equations (55a) and (55b). A detailed example of the different boundary condition implementations in FEniCS is shown in Appendix A.
For the wave equation in PH form, setting p at the boundary is a Dirichlet condition, equivalent to setting ρ w e p , and setting (q · n) is a Neumann condition, equivalent to setting 1 kw (e q · n). The reason that we can implement Dirichlet and Neumann conditions in the weak form is because we integrate both Equations (13a) and (13b) by parts, giving boundary terms for Dirichlet and Neumann conditions in Equations (55a) and (55b). This approach differs from the PFEM of Cardoso-Ribeiro [10], where integration by parts is only used on a subset of the governing equations. Our approach thus has the advantage of allowing mixed boundary conditions without the need for Lagrange multipliers. Therefore, our method results in a matrix system of equations that can be solved as an ordinary differential equation (ODE), rather than a differential algebraic equation (DAE), which is easier to solve in general. All methods in this section use a time step of t = 5 × 10 −4 . Figure 3 shows the resulting Hamiltonian,Ĥ , for the input in Equations (52a)-(52c). The Hamiltonian is expected to be constant after 0.25 s, because the input boundary condition is set to zero, therefore, no energy flows into or out of the domain.
In Fig. 3(a) we show that the implicit Euler (IE) scheme incorrectly dissipates energy and the explicit Euler (EE) scheme has a fictitious increase in energy, resulting in instability. At t ≈ 0.9 s the SE integration schemes for strong and weak boundary conditions show an undesirable non-constant Hamiltonian. A zoomed-in image of the non-constant behaviour is also shown in Fig. 3(b). The bump at t ≈ 0.9 s is likely due to there being high-order dynamics when the wavefront approaches the right boundary that are not accounted for in the first-order SE scheme. To remedy this situation, three second-order integration schemes are implemented, Explicit Heun's (EH) (also called improved Euler [39]), SV and SM. All of these methods drastically decrease the bump at t ≈ 0.9 s. It should be noted that decreasing the time step of the first-order methods also has the same effect of decreasing the bump (not shown). Interestingly, EH removes the bump completely, however, SV does not. SM also removes the bump completely because the SM scheme conserves quadratic invariants exactly for linear systems [20]. Figure 3(b) shows that for SE and SV the Hamiltonian oscillates about a conserved average Hamiltonian. This plot also shows that SV has an approximately constant Hamiltonian (apart from the aforementioned bump at t ≈ 0.9 s), whereas EH has a gradually increasing Hamiltonian, indicating that energy is not conserved for long times. For SV, however, the Hamiltonian oscillates about a conserved value and therefore, is conserved for long times. The bound on the oscillations also converges towards the conserved value as the time step is decreased [20]. The energy bound for SE and SV is likely higher at the bump due to the reflecting boundary. This would hint that the non-perfect energy conservation of SE and SV is heavily influenced by the boundary and not so much by the internal domain. The SM One problem of the first-order SE weak method in Fig. 3 is its oscillatory behaviour. These oscillations may be due to applying an essential boundary condition (EBC) weakly, without a penalty method. Typically EBC's are applied with a penalty method such as Nitsche [25,30], however, according to Scovazzi, penalty methods are not required due to the hyperbolic nature of the wave equation [37]. This could mean that the oscillations may simply be the expected oscillations of low order symplectic time integration schemes. As anticipated by the proof of boundedness in [18], the oscillations are reduced when the time step is decreased (not shown) or the order of the symplectic method is increased, as seen with the SV method.
To compare the energy conservation of weak and strong boundary condition implementations, we define the energy residual at time t f as whereĤ , which is equal to the internal energy in the domain at time t f , is calculated from Equation (27) withq andp variables at time t f . Kotyczka et al. [23] shows that implicit Gauss-Legendre schemes such as SM conserve the discrete energy exactly for linear PH systems. However, for schemes such as SV and SE that do not conserve energy perfectly, we expect an energy error from both the non-conservativity within the domain and a nonexact energy transferred through the boundary. Therefore, the second term in Equation (59) is simply the energy that we expect to have transferred through the left boundary out of the domain (the term is negative for flow into the domain) between the initial time t 0 and time t f . Since we want to calculate the numerical energy flow out of the boundary, care must be taken to evaluate the output energy calculation in the same way as the numerical time integration scheme. This ensures that there is no discrepancy between the accuracy in which the two terms in Equation (59) are calculated. To ensure this consistent numerical accuracy, variables with a overscript are the effort or state variables evaluated at time steps that correspond with the chosen time integration scheme. For example, the expected energy contribution of the SV scheme at the boundary can be calculated by decomposing the SV method, as shown in a concise form here, into two adjoint SE method steps with half timestep each [20], this gives where F () denotes a function of the entries in the bracket. Then, by knowing that the superscripts of q and p L are the same as the superscripts of q and p in Equation (61a) for the first half time step, and the same as in Equation (61b) for the second half time step. We can calculate the residual in two half time steps, as follows from Equation (59), The energy residual in Equation (59) is plotted for multiple time integration schemes in Fig. 4.
It is clear from the bounded energy residual of the SE and SV schemes and the exactly conserved energy of the SM scheme in Fig. 4 that the weak boundary condition implementation conserves energy for long times. This agrees with our proof of Theorem 2.2.1 and the expected bounded energy residual of symplectic time integrations schemes. Again, for the SE and SV schemes, the energy is conserved in an average sense and oscillations about the conserved energy do occur. As expected, the SM scheme conserves energy exactly, with an energy residual of < 10 −12 , which is round-off error. The strong implementation of the input Dirichlet boundary condition does not conserve energy. More exactly, the strong implementation's energy residual is dependent on the refinement of the spatial mesh. To display this effect, the energy residual of the wave equation with SV time integration is analysed for varying element characteristic length ( √ Mean Element Area) and is shown in Fig. 5. A quadratic trend is plotted to show that the energy residual when using a strong boundary condition implementation has a quadratic dependence on the element characteristic length. This differs from the weak boundary condition implementation, which conserves energy independently of the mesh refinement. This agrees with our proof of Theorem 2.2.1, i.e., that our spatial discretisation is perfectly energy-conserving. Since the energy error is bounded for symplectic time integration schemes [19], the energy error for the weak boundary implementation is only dependent on the step size of the symplectic time integration scheme.
Temporal integration schemes cannot, in general, conserve both the exact energy and the symplectic structure of the system [20]. However, a general result that applies to nonlinear systems, as well as our linear system, is that conserving the symplectic structure results in a bounded energy error which decreases as the time step is reduced. In subsequent chapters, we focus on the results of the SV scheme rather than the SM scheme to show energy conservation results that resemble what we expect for general non-linear problems.

Wave Equation Comparison with Analytical Solution
In this section, our numerical model with a weak boundary condition implementation is compared against an analytical solution to ensure stable spatial and temporal convergence  [32]. The initial and boundary conditions are given as q(x, L y , t) · n = 0 .
The analytical solution is then given by where all constants are given in Appendix D.

Spatial Convergence
The error between the numerical model and the analytic solution for a range of characteristic element lengths is evaluated to show the spatial convergence of the model. A table of the number of elements, with corresponding error and convergence details for each refinement level, is shown in Appendix E. The SV time integration scheme was used with a time step of 5 × 10 −4 s. The L 2 error norm for each step n is defined as where p a (t n ) is the exact solution evaluated at the step n. This error is calculated accurately with the 'errornorm' function in FEniCS. The maximum L 2 error norm over 1.5 s of simulation is plotted in Fig. 6 for each characteristic element length. In the figure legend, the numbers following P and RT denote the order of the Lagrange and Raviart-Thomas elements, respectively, i.e., P1RT2 uses first order Lagrange elements and second order Raviart-Thomas elements. As can be seen, the L 2 error norm for all element order combinations shows standard convergence against characteristic element length. Here we define standard spatial convergence as convergence of order O(( x c ) k+1 ), where k is the order of the methods lowest order basis function and x c is the characteristic element length.

Modal Analysis
To ensure the correct handling of mixed boundary conditions, the eigenvalues of the model were verified to be accurate by performing a modal analysis. The analytic eigenvalues were calculated by separation of variables of Equation (1) into This method can simply be shown to give the following three first order eigenvalue problems, where λ and μ are eigenvalues of the first order spatial systems and ω 2 = c 2 (μ + λ) is the eigenvalue or the squared eigenfrequency that we want to predict. For the boundary conditions in Equations (63a)-(63e), the real components of the eigenvalues in Equations (67b) and (67c) are zero and the complex components are Therefore, from Equation (67a), the analytic eigenfrequencies we wish to find are given by To predict the eigenfrequencies of our model, Equations (55a) and (55b) were discretised in FEniCS, creating the following system of equations, where, M e , K e , and L e are the mass, stiffness, and boundary matrices which determine the eigenfrequencies of the system. The M −1 e (K e + L e ) matrix was then input into NumPy's eigensolver to calculate the eigenfrequencies of the system. A plot of the complex part of the first 50 analytic and modelled eigenfrequency pairs for a mesh with 1322 elements is shown in Fig. 7. The real parts of all eigenfrequencies equal zero, as expected for the wave equation with no damping.
As shown in Fig. 7, the eigenfrequencies are predicted accurately with no spurious modes. A plot of the convergence of the eigenfrequencies with respect to the characteristic element length is shown in Fig. 8. This shows standard quadratic convergence for first order Lagrange and Raviart-Thomas elements.

Time Convergence
The numerical error for this model is heavily dominated by the spatial discretisation. Therefore, to show time convergence, third order Lagrange and Raviart-Thomas elements are used. This ensures that the initial spatial error is small, therefore, the error that propagates through the domain is due to the temporal discretisation. A convergence study is done for a 1.5 s simulation and an error growth study is assessed for long-times (t = 10 s). The long-time analysis shows the rate at which the error grows through time and displays our models effectiveness for long-time simulations. Although the rate of error growth for different symplectic schemes is well known [18], this section importantly shows that our spatial discretisation does not deteriorate the expected error growth rate.
The SE and SV time integration schemes have been implemented to show the temporal convergence of our method. Linear and quadratic convergence is shown for the SE and SV schemes, respectively in Fig. 9, where the maximum L 2 error over the 1.5 s simulation is plotted for a range of time steps. This shows that our method gives standard temporal convergence, where we define standard temporal convergence as convergence of order O(( t) r ), where r is the order of the time integration scheme. These simulations were performed on a spatial discretisation with 9358 elements.
The order at which the state variable error of a method grows is a common metric for the accuracy of symplectic and multi-symplectic methods, as assessed extensively in Hairer's book [20]. Hairer shows that symplectic methods have state variable error growth of order O(t ( t) r ), where r is the order of the time integration scheme. The time step convergence in Fig. 9 for fixed final time t = 1.5 s confirms that the error converges with order O(( t) r ) Therefore, observing an error growth proportional to t when simulated for long times is sufficient to show the correct order of error growth, O(t ( t) r ). The L 2 error norm for a 10 s simulation of the wave equation is shown in Fig. 10 for various time integration schemes. To give a fair comparison, the schemes have time steps that result in the same number of function evaluations, 0.00025 s for IE and SM and 0.0005 s for SV. To display a result that does not blow up immediately, the time step for EH is decreased even further to 0.000125 s.
In Fig. 10 the EH method blows up with exponential error growth because it does not have a bounded energy residual. This behaviour is typical of fully explicit methods, which can be unstable for long times [19]. The IE method has large error growth due to its inherent energy dissipation. This error increase tapers off (not shown) due to the complete loss of energy in the numerical model, this results in p n approaching a constant zero throughout the domain. The symplectic methods both show an error that is linearly dependent on time, as required to validate that the error growth is proportional to O(t ( t) r ). The high-frequency oscillations of the L 2 error norm in Fig. 10 are due to p n varying from being zero throughout the domain to having a maximum, as shown in Fig. 11. Therefore, the oscillation frequency is the frequency that p oscillates from maximum/minimum to zero.

Electromechanical Lumped Parameter Model
This section details a simple LPM for a linear actuated electric motor. The system diagram and bond-graph schematic are shown in Fig. 12, with constants defined in Appendix D. The current of the electrical system is i and the displacement of the linear motor is s. The bondgraph methodology is a modular approach for LPMs that ensure conservation of energy within and between models, for a review see the work by Gawthrop [16]. The PH framework extends from bond-graphs to also allow continuum models that conserve energy.
The Hamiltonian for this system is where the canonical momentum for the mechanical subsystem is given by h M = mṡ and the electrical system equivalent of the canonical momentum is the magnetic flux linkage, Fig. 12 (a) System diagram and (b) bond-graph diagram for an electromechanical system, where Se, R, I, C, and GY denote effort sources, dissipative components, inductive/mass storage components, capacitive/spring storage components, and gyrator components respectively denoted by h E = L E i. The canonical PH form for this system is where B u is the input control matrix with corresponding input, u and output, y u . Also, B F is the boundary port matrix with corresponding boundary force F b and boundary output y F . Including the constitutive laws and evaluating the B u and B F matrices gives In this paper, the resistances, R E and R M , are set to zero, as a non-dissipative system is required to display a conserved Hamiltonian. The values of all constants used in the simulations are given in Appendix D. From Sect. 6 onwards F b will be the reaction force from the coupled wave equation. The first part of Equation (73) can be written as a typical linear system of ODE's with a state vector y = [h E , h M , s] T and a control/interconnection vector u = [u, F b , 0] T , asẏ To create a monolithic coupling of this LPM with the wave equation from Sect. 3 the LPM needs to be implemented in FEniCS. The ODE in Equation (74) can be implemented in FEniCS by using the 'real function space', which assumes a function has one value over a domain, i.e., it has no spatial dependence. This makes the domain that the equations are implemented on irrelevant to the ODE. Equation (74) can then be implemented by multiplying by a trial function, v y , and taking the trace against the border of an arbitrary domain. To allow easier coupling in the following section the trace is taken against the left boundary of the wave domain. The remainder of this section details the implementation of the LPM in FEniCS. Although in the following sections we use an SV time integration scheme, the SE scheme is detailed here to provide easier understanding to the reader. For the implementation of the SV method see the GitHub repo. https://github.com/FinbarArgus/portHamiltonian_FEM.git For the SE scheme, the ODE can be implemented by solving the matrix systems of equations generated by at each time step. This method of implementing a LPM by assigning the variables as real functions over the whole domain is not optimally efficient and thus, future work should look at developing a method specifically for LPMs that is compatible with FEniCS. y m is a vector of state variables at the previous time step and y a is a vector of state variables at a combination of current and previous time steps, as shown here, Once real function spaces are created and the vectors in Equation (76) are formed, the variational form in Equation (75) can be expressed and solved with the Python code in Appendix B.

Coupling with the Electromechanical Model
In this section, the wave equation from Sect. 3 is coupled with the EM PH model from Sect. 5. A schematic of the combined model is shown in Fig. 13 for a rectangle wave domain.
Following the previous sections, the boundaries denoted ∂ M and ∂ R have zero Neumann and zero Dirichlet conditions, respectively. The ∂ L boundary is the Dirichlet boundary connection with the LPM model. All boundary conditions are implemented in a weak manner. Since the only boundary connection is a Dirichlet boundary, the canonical form of Equation (49) can be used to set up the interconnection with Equation (73). The relationship between inputs and outputs of the wave domain and EM domain can be written as where W is a compact operator and W * is its adjoint operator. The duality pairings for the inputs and outputs are where the ·|· ∂ L duality product is an L 2 inner product that acts over the connection boundary and the ·, · inner product is the standard inner product in R. We know that the velocity of the left wave boundary is directly set by the output velocity of the EM system, therefore, where w transforms a scalar into a constant function over the boundary, with the value of y F . Substituting the boundary output relation from Equation (73) gives In practice, to turn the scalar value into a function it must be turned into a constant vector and then the vector must be dot producted with a vector of basis functions. Knowing that u p = θ T pû p we get, Now, we enforce that every point on the left boundary of the wave domain has the same vertical velocity (ê p | ∂ L ) as the vertical velocity of the output motor ( h M m ). Therefore,û p = −ê p | ∂ L has the value of − h M m for each of its DOFs, which gives The following energy conserving relation between the L 2 and R inner products in Equation (78) can be used to determine W * , Evaluating the left-hand side of Equation (84) gives Equating to the right-hand side of Equation (84) gives which results in Substituting y q from Equation (35) and e q = k w q gives Therefore, Equations (80) and (88) are the energy conserving interconnection relations between the two domains. The total Hamiltonian of the combined system is given by Finally, the canonical forms of Equations (47) and (73) can be combined to give the input-state-output PH form of the interconnected system, ⎡ where D p and D 0 p are the interconnection matrices of the skew symmetric system matrix. For skew-symmetry of the system matrix to hold the following must be true To verify that this skew symmetry holds we calculate D p and D 0 p by discretising Equations (80) and (88). Firstly, Equation (80) is discretised by multiplying on the left by the trial basis function, ψ q , and integrating over the boundary, as was done in the formulation of the boundary term in Equation (19), where the final step comes from the definition of B pb in Equation (48). Equating this with the boundary term in the second row of Equation (90) gives To discretise Equation (88), we first substitute Equation (83) and apply w = w T , because w is a constant function, Discretising q · n in the same way as e q · n was discretised in Equation (18) gives, where 1 T simply adds up the force contribution from each boundary DOF to calculate F b . Equating the boundary term in the fourth row of Equation (90) with the boundary term in the second row of Equation (73) gives Therefore, the discrete interconnected system retains a skew symmetric matrix, as required for a port Hamiltonian system. Equation (90) encompasses the canonical form, the dynamic equations and the constitutive law equations of the interconnected system.
It is important to note that if the right-side boundary condition, p R , was non-zero then there would be an extra boundary input and output in the canonical form of Equation (90). Also, if the top and bottom Neumann boundaries were non-zero conditions, the canonical form of Equation (90) would have interconnection inputs u int p and outputs y int q with the same definitions as the discrete version of Equation (35). A canonical form for the Neumann conditions would also have to be created with inputs/outputs for the boundary u q , y p and for the interconnection u int q , y int p with the same definitions as the discrete version of Equation (34). The two canonical forms could then be connected with Equation (50). Although the formulation of the interconnected canonical form seems complicated, it is not necessary for implementation. The full system canonical form is only discussed to reassure the reader that the total system is equivalent to a combination of canonical forms and therefore is a Stokes-Dirac structure and can have mixed causality boundary connections. The usefulness of the canonical forms, Equations (47) and (73), is that they allow calculation of the interconnection relations, which can then be implemented in FEniCS, as shown in Appendix C.

Interconnection Model Results
This section displays the results for the coupled wave-EM model for a sinusoidal input voltage. The time step for the simulations is 5×10 −4 s, using an SV time integration scheme. The problem is solved on both a simple rectangle domain with 21110 elements and a square domain with central input and 13067 elements.

Rectangle Domain
The input voltage condition for the rectangle domain is given by The p variable of the wave equation is displayed after 0.3 s and 1.1 s in Figs. 14 and 15, respectively. The initial large sinusoidal wave, which is due to the input over the first 0.25 s, flows through the domain as expected. This wave is followed by smaller repeating waves caused by the lingering oscillations of the EM system.
To confirm that the energy error of the model is bounded, the energy residual of each domain and the total energy residual is plotted in Fig. 16. The residual for each domain is Wave residual is the difference between the energy in the wave domain and the accumulated energy that has entered the wave domain through its boundaries. EM residual is the difference between the energy in the EM domain and the accumulated energy that has entered the EM domain through its boundaries the difference between the accumulated energy that has entered the domain and the internal energy in the domain.
As shown, the energy residual is bounded for increasing times, as expected for a symplectic time integration scheme. For a second-order symplectic method such as SV, the energy residual bound should be quadratically dependent on the time step, this relationship is proven in [18] and expressed more generally as To ensure that the bounded energy residual is indeed quadratically dependent on the time step, the energy residual maximum over a 20 s simulation is plotted against the step size in Fig. 17. As discussed in Sect. 4 The SM scheme can also be used to conserve quadratic invariants exactly. As shown in Fig. 17 we achieve an energy residual of < 10 −11 for a 20 s simulation with SM, confirming that our model, with an SM scheme, conserves energy exactly for linear, coupled LPM-continuum systems.

Square Domain with Central Input Boundary
A diagram for the domain of this section is shown in Fig. 18 and the input voltage is the same as Equation (98). Following the notation of the previous sections, the boundaries denoted At long times the uniform circular wave structure breaks down due to waves rebounding off the walls. However, the wave behaviour should still retain some structure, for example, there should be symmetry about the midline parallel to the x axis. Figure 21 shows p n after 8 s of simulation, once the uniform wavefronts have completely broken down. As can be seen, there is still symmetry about the midline parallel to the x axis, further showing this methods ability to accurately model the physical structure of the system. The energy residual for an 8 s simulation with the SV scheme is shown in Fig. 22. As in previous sections, the energy residual is oscillatory and bounded, as expected. Also, as in the previous sections, when using the SM method the energy residual is < 10 −12 , therefore, exactly conserved to within round-off error.

Conclusion
For the modelling of continuum systems, a port-Hamiltonian, Galerkin finite element method has been developed that, in general, has a bounded energy residual and linear longtime error growth for the state variables. This method allows mixed boundary conditions without the need for Lagrange multipliers or user-defined parameters. The discretisation is shown to be equivalent to a coupling of canonical port-Hamiltonian forms that allows Fig. 22 Energy residual vs time for the interconnected wave-EM simulation with a square domain, using Störmer-Verlet time integration. Total Residual is the sum of the energy residuals from the wave and EM domains. Wave residual is the difference between the energy in the wave domain and the accumulated energy that has entered the wave domain through its boundaries. EM residual is the difference between the energy in the EM domain and the accumulated energy that has entered the EM domain through its boundaries mixed interconnections with other canonical port-Hamiltonian models. The discretisation is also shown to be symplectic in both time and space. For our specific 2D linear wave equation system we also show exact energy conservation with a symplectic (implicit) midpoint method that guarantees conservation of quadratic first integrals. We also compare against an analytical solution and show standard order of convergence for the state variables with respect to the temporal and spatial discretisation. A modal analysis is performed and the eigenvalues are verified to be accurate. The boundary conditions are implemented in variational form, in a weak manner, without the need for penalty methods. In addition to this, the method is capable of monolithic coupling with arbitrary LPMs. The coupled model is shown to also have a bounded energy residual with a standard temporal order of convergence for the SV time integration scheme and exact energy conservation for the SM time integration scheme. The example model of a 2D linear wave equation coupled with an EM linear actuator is a good proof of concept for more advanced couplings between Hamiltonian PDEs and LPMs. Future work will be done on implementing control algorithms for coupled models, in order to improve control of multiphysics, multidomain, and non-linear problems.

Appendix A: FEniCS Code Example
To give the reader details on how the variational forms in Sect. 3 are implemented in FEn-iCS, this section shows code snippets for the SE scheme with weak form and strong form boundary conditions. The full code can be found at https://github.com/FinbarArgus/portHamiltonian_FEM.git Firstly FEniCS and the meshing software mshr are imported, then the domain shown in Fig. 2 is created with the following lines of code, from f e n i c s i m p o r t * i m p o r t mshr m a i n R e c t a n g l e = mshr . R e c t a n g l e ( P o i n t ( 0 . 0 , 0 . 0 ) , P o i n t ( L_x , L_y ) ) mesh = mshr . g e n e r a t e _ m e s h ( m a i n R e c t a n g l e , r e s ) where res is the resolution of the mesh. The function spaces, trial functions, test functions and functions for the DOFs at the previous and current time steps can then be created as follows, (100) The variational form in Equation (75) that allows the EM domain input boundary term, Equation (87), to be applied in FEniCS is where v y2 is the 2nd component of v y , the ODE test function vector. Substituting Equation (87) into Equation (101) and taking v y2 out of the integral, since it is a real function space that is constant over the boundary, gives Now, since the inside integral evaluates to a constant over the boundary, the outside integral can be evaluated to give v y2 L y ∂ L y q ds L .
Finally, substituting y q from Equation (35) and e q = k w q gives k w L y v y2 ∂ L (q · n)ds L .
The interconnection with boundary interconnection forms, Equations (100) and (104), can be implemented in FEniCS for the SE scheme with the following code,  Wave domain stiffness k w 3.0 k g m s −2 Wave speed squared c 2 (k w /ρ w ) 1.5 m 2 s −2 Wave domain x-length