Hassan Fahs
email: hassan.fahs@gmail.com

High-order discontinuous Galerkin method for time-domain electromagnetics on non-conforming hybrid meshes

Keywords: Maxwell's equations, discontinuous Galerkin method, hybrid meshes, non-conformal meshes

We present a high-order discontinuous Galerkin (DG) method for solving the time-dependent Maxwell equations on non-conforming hybrid meshes. The hybrid mesh combines unstructured tetrahedra for the discretization of irregularly shaped objects with a hexahedral mesh for the rest of the computational domain. The transition between tetrahedra and hexahedra is completely non-conform, that is, no pyramidal or prismatic elements are introduced to link these elements. Within each mesh element, the electromagnetic field components are approximated by a arbitrary order nodal polynomial and a centered approximation is used for the evaluation of numerical fluxes at inter-element boundaries. The time integration of the associated semi-discrete equations is achieved by a fourth-order leap-frog scheme. The method is described and discussed, including algorithm formulation, stability, and practical implementation issues such as the hybrid mesh generation and the computation of flux matrices with cubature rules. We illustrate the performance of the proposed method on several two-and three-dimensional examples involving comparisons with DG methods on single element-type meshes. The results show that the use of non-conforming hybrid meshes in DG methods allows for a notable reduction in computing time without sacrificing accuracy.

Introduction

Nowadays, a variety of modeling strategies exist for the computer simulation of electromagnetic (EM) wave propagations in the time-domain. In particular, the ability to deal with complex geometries has been extensively studied in the context of unstructured methods as alternatives to the well established FDTD method [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF]. In spite of its versatility, the conventional FDTD method suffers from serious accuracy degradation when modeling curved objects and small geometrical features. This is due to the commonly used staircase approximation (Cartesian grid) which introduces large discretization errors, even when the grid size is very small. In contrast, unstructured solvers such as the finite element time domain (FETD) [START_REF] Pernet | High spatial order finite element method to solve Maxwell's equations in time domain[END_REF] or the finite volume time domain (FVTD) [START_REF] Piperno | A nondiffusive finite volume scheme for the three-dimensional Maxwell's equations on unstructured meshes[END_REF] methods can easily handle complex boundaries by dealing with unstructured tetrahedral meshes. Both methods are more expensive than the FDTD method in terms of computing time and memory requirement. This motivates the development of hybrid solution strategies, in which the computational efficiency of the Cartesian grid is combined with the geometric flexibility of the unstructured mesh. The FDTD is then used in large homogeneous volumes of the computational domain while the FVTD or FETD methods are employed in small regions near small scale features. For instance, hybrid FETD-FDTD methods have been developed [START_REF] Abenius | Hybrid time domain solvers for the maxwell equations in 2d[END_REF][START_REF] Wu | Hybrid finite-difference time-domain modeling of curved surfaces using tetrahedral edge elements[END_REF], but they suffer from spurious reflections and late-time instabilities, which limit their usefulness. A stable hybrid FETD-FDTD method has been described in [START_REF] Rylander | Stability of explicit-implicit hybrid time-stepping schemes for Maxwell's equations[END_REF]. To achieve the time-stability, the scheme requires an implicit time integration in the unstructured elements but also a transition layer of pyramids between the FDTD cubes and tetrahedra to ensure the mesh conformity. Low order hybrid FDTD-FVTD methods are proposed in [START_REF] Edelvik | Explicit hybrid time domain solver for the Maxwell equations in 3D[END_REF][START_REF] Ferrieres | Application of a hybrid finite difference/finite volume method to solve an automotive EMC problem[END_REF]. These methods require the use of an intermediate layer of hexahedra between the Cartesian grid and the tetrahedral mesh. Each cell in the layer is treated by the FVTD method, which offers a correspondence between a quadrilateral face and two triangles in terms of the exchange of fluxes. Again, weak instabilities have been reported for long-time simulations.

The last few years have seen an increasing interest in so-called discontinuous Galerkin time domain (DGTD) methods. The DG methods are a class of FE methods based on completely discontinuous piecewise polynomial spaces for the numerical solution and the test functions. To obtain highly accurate and stable DG methods, suitable numerical fluxes need to be designed over elemental interfaces. We refer to the lecture notes [START_REF]Discontinuous Galerkin Methods[END_REF] and the textbook [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: Algorithms, Analysis and Applications[END_REF] for details and history about DG methods. In particular, DG methods can easily handle elements of various types and shapes [START_REF] Gassner | Polymorphic nodal elements and their application in discontinuous Galerkin methods[END_REF], irregular non-conforming meshes [START_REF] Dosopoulos | Non-conformal and parallel discontinuous galerkin time domain method for maxwell's equations: EM analysis of IC packages[END_REF][START_REF] Fahs | Preliminary investigation of a non-conforming discontinuous Galerkin method for solving the time-domain Maxwell equations[END_REF], and even locally varying polynomial degree [START_REF] Fahs | Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic wave propagation[END_REF], and hence offer great flexibility in the mesh design. They also lead to (block-) diagonal mass matrices and therefore yield fully explicit, inherently parallel methods when coupled with explicit time schemes. Furthermore, DG methods are flexible with regards to the choice of the time-stepping scheme [START_REF] Fahs | Investigation on polynomial integrators for time-domain electromagnetics using a high-order discontinuous Galerkin method[END_REF]. The DG methods have now acquired a level of maturity and have successfully penetrated several scientific and technological communities following to their adaptation to increasingly complex modeling contexts [START_REF] Busch | Discontinuous Galerkin methods in nanophotonics[END_REF][START_REF] Fahs | Calculation of the SAR induced in head tissues using a high order DGTD method and triangulated geometrical models[END_REF][START_REF] Schnepp | Efficient large scale electromagnetic simulations using dynamically adapted meshes with the discontinuous Galerkin method[END_REF]. It is worthwhile to note that the DG method has also been adopted for the first time in a commercial software as the time domain alternative of a very well known EM wave simulation tool [START_REF] Songoro | Keeping time with Maxwell's equations[END_REF]. While DG methods on tetrahedral meshes have been extensively studied [START_REF] Dosopoulos | Non-conformal and parallel discontinuous galerkin time domain method for maxwell's equations: EM analysis of IC packages[END_REF][START_REF] Fahs | Calculation of the SAR induced in head tissues using a high order DGTD method and triangulated geometrical models[END_REF][START_REF] Hesthaven | Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell's equations[END_REF], their development on hexahedral meshes has received relatively less attention [START_REF] Cohen | A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell's equations in time domain[END_REF]. One reason is that tetrahedral meshing is highly automated while high quality hexahedral meshing of complex geometries commonly requires user intervention and is labor intensive. This motivates the development of DG methods on conforming hybrid meshes consisting of tetrahedra and hexahedra with transition layers of prisms or pyramids [START_REF] Bergot | Higher-order discontinuous Galerkin method for pyramidal elements using orthogonal bases[END_REF][START_REF] Sevilla | The use of hybrid meshes to improve the efficiency of a discontinuous Galerkin method for the solution of maxwell's equations[END_REF]. Despite these advances, the exploitation of non-conforming hybrid meshes in the context of a DG-based numerical methodology applied to wave propagations, has been the subject of a very few studies so far. A 2D contribution in this context is the work of Durochat et al. [START_REF] Durochat | DGTD method on hybrid meshes for time domain electromagnetics[END_REF][START_REF] Durochat | High order non-conforming multi-element discontinuous Galerkin method for time-domain electromagnetics[END_REF] who proposed a DG method formulated on hybrid rectangular-triangular meshes for Maxwell's equations.

We are concerned here with the design of a high-order nodal DGTD method for solving Maxwell's equations on non-conforming hybrid meshes in 2D and 3D. The proposed method is based on centered numerical fluxes and a fourth-order leap-frog time integration scheme. The hybrid mesh adopted here combines tetrahedral elements for the discretization of irregularly shaped objects with hexahedral elements to fill the rest of the domain. Thus, the resulting mesh is non-conformal in the sense that hanging nodes are allowed and no additional element types (e.g., pyramids or prisms) are needed to join tetrahedra with hexahedra. The non-conformal connection between these elements is ensured thanks to the correct computation of flux terms by using efficient cubature rules over polygons. The motivation behind this work stems partly from the following observations: first, the final mesh of an object is made of separately constructed meshes and a conforming assembling is a very hard task. A good example to illustrate this point is the interior penalty DGTD method formulated on non-conforming tetrahedral meshes presented in [START_REF] Dosopoulos | Non-conformal and parallel discontinuous galerkin time domain method for maxwell's equations: EM analysis of IC packages[END_REF] for the EM simulation of 3D integrated circuit packaging. For the complicated geometries involved in such applications, the construction of a conforming mesh is a very hard task. Instead, a subdomain-based meshing approach is particularly attractive and allowing meshes with different element-types is highly desirable in this context. Second, regarding the computational effort, it turns out to be more efficient to compute on hexahedral meshes instead of tetrahedral meshes to reach a desired error level. Therefore, the performance of the DG method can be increased by combining tetrahedra with hexahedra. The outline of the paper is as follows. In Section 2 we briefly recall the basic features of DGTD discretizations based on totally centered numerical fluxes and formulated on single element-type meshes composed of triangles or quadrilaterals in 2D and tetrahedra or hexahedra in 3D. Then, we formulate the DG method on nonconforming hybrid triangular-quadrilateral or tetrahedral-hexahedral meshes in Section 3. Subsequently, we analyze the stability of the resulting algorithm and establish a sufficient CFL-like condition. Next, in Section 4, we discuss some practical implementation issues including the generation of non-conforming hybrid meshes and the computation of the flux matrices at hybrid interfaces between adjacent elements of different types. Finally, we present numerical results for several test-cases in Section 5, and draw some conclusions of the study in Section 6.

Governing equations and DGTD formulation

We consider the time-dependent Maxwell equations for linear, isotropic, non-dispersive media. The electric field E and the magnetic field H verify

ε(x)∂ t E(x, t) = ∇ × H(x, t), µ(x)∂ t H(x, t) = -∇ × E(x, t), (1)
where ε and µ are respectively the permittivity and the permeability of the medium. These equations are posed on a bounded domain Ω of R d , d ≥ 2. At the domain boundary ∂Ω, we apply either perfect electric conductor (PEC) conditions (n × E = 0) or absorbing boundary conditions (n

× E + cµn × (n × H) = 0).
Here n denotes the unit outward normal to ∂Ω and c = (εµ) -1/2 is the speed of propagation.

We assume we dispose of a regular or irregular partition T h of Ω into n e non-overlapping elements κ i , which can be triangles or quadrilaterals in 2D and tetrahedra or hexahedra in 3D. Thus a (d -1)-dimensional face of each element κ i in T h is allowed to contain hanging nodes. We shall suppose that each κ i ∈ T h is the image of a fixed master element κ r , under a smooth bijective mapping

Ψ κ i : κ r ∋ ξ ξ ξ → x ∈ κ i , with nonsingular Jacobian J κ i = ∂x/∂ξ ξ ξ, that is, κ i = Ψ κ i (κ r ) for all κ i ∈ T h , where κ r is either the open unit simplex T d = {ξ ξ ξ = (ξ 1 , . . . , ξ d ) ∈ R d : ξ i ≥ 0, d i=1 ξ i ≤ 1} or the open unit hypercube D d = {ξ ξ ξ = (ξ 1 , . . . , ξ d ) ∈ R d : 0 ≤ ξ i ≤ 1} in R d .
Within this construction we allow for elements κ i with curved edges or faces. On the master element κ r , with ξ ξ ξ = (ξ 1 , . . . , ξ d ) ∈ κ r and α α α = (α 1 , . . . , α d ) ∈ N d 0 , we define spaces of polynomials of total degree ℓ ≥ 1 as follows

Q ℓ (D d ) = span{ξ ξ ξ α α α : 0 ≤ α i ≤ ℓ, 1 ≤ i ≤ d}, P ℓ (T d ) = span{ξ ξ ξ α α α : 0 ≤ |α α α| ≤ ℓ}.
The dimension of each of these spaces depends on the spatial dimension d and on the degree ℓ and is given by

N(ℓ, d) =        (ℓ + 1) d if κ r = D d , (ℓ + d)!/ℓ!d! if κ r = T d .
As a particular basis for Q ℓ (D d ) and P ℓ (T d ) consider the fundamental polynomials {L (ℓ) k } N ℓ k=1 of multivariate Lagrange interpolation, corresponding to a set of nodal points S ℓ = {η j , j = 1, . . . , N ℓ }, such that L (ℓ) k (η j ) = δ jk for η j ∈ S ℓ . The distribution of the nodes η j is a key issue for the properties of the interpolation, especially for very high-order approximations. For the unit simplex T d , we choose the nodal set derived from the electrostatics principle [START_REF] Hesthaven | From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex[END_REF][START_REF] Hesthaven | Stable spectral methods on tetrahedral elements[END_REF], while the Gauss-Legendre-Lobatto points are used for the hypercube D d .

To each κ i ∈ T h we assign an integer ℓ i ≥ 1 that is the local polynomial degree inside κ i ; collecting the ℓ i and Ψ κ i in the vectors ℓ ℓ ℓ = {ℓ i : κ i ∈ T h } and Ψ Ψ Ψ = {Ψ κ i : κ i ∈ T h }, respectively, we introduce the finite element space

V p (Ω, T h , Ψ Ψ Ψ) = u ∈ [L 2 (Ω)] d : u |κ i • Ψ κ i ∈ [Q ℓ i (κ r )] d if Ψ -1 κ i (κ i ) = D d ,
and

u |κ i • Ψ κ i ∈ [P ℓ i (κ r )] d if Ψ -1 κ i (κ i ) = T d , ∀κ i ∈ T h ,
where [L 2 (Ω)] d is the space of square-integrable functions on Ω. For two distinct elements κ i and κ k in T h , the intersection a ik = κ i ∩ κ k is a (d -1)-dimensional face in T h , with unitary normal vector n ik , oriented from κ i towards κ k . For the boundary faces, the index k corresponds to a fictitious element outside the domain. Finally, we denote by F I and F B the union of all interior and boundary faces of T h , respectively.

The derivation of the DG scheme for the simulation of electromagnetic wave propagation on purely simplicial meshes has already been published in previous works [START_REF] Fahs | Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic wave propagation[END_REF][START_REF] Fahs | A high-order non-conforming discontinuous Galerkin method for time-domain electromagnetics[END_REF]. Therefore, we avoid the repetition of a comprehensive derivation and refer the reader to these works for details. In brief, the electric and magnetic fields are approximated inside each finite element κ i by a linear combination of basis functions ϕ i j = L (ℓ i ) j • Ψ -1 κ i of degree ℓ i with support κ i and with time-dependent coefficient functions E i j (t) and H i j (t) as follows

E i (x, t) = 1≤ j≤N i E i j (t)ϕ i j (x), H i (x, t) = 1≤ j≤N i H i j (t)ϕ i j (x). (2)
Here, the index j indicates the jth basis function and N i = N(ℓ i , d) is the number of degrees of freedom inside κ i . We denote by Ēi and Hi , respectively, the column vectors (E i j ) 1≤ j≤N i and (H i j ) 1≤ j≤N i . As usual for DG schemes, the Maxwell system, Eq. ( 1), is multiplied by a test function ϕ ∈ span{ϕ i j , 1 ≤ j ≤ N i } and integrated over each κ i . After integration by parts, inserting the DG approximation, Eq. ( 2), and after applying a centered numerical flux in the boundary integrals, the semi-discrete formulation of the scheme in the physical element κ i reads as

M ε i ∂ t Ēi = K i Hi - a ik ∈F I S ik Hk - a ik ∈F B S ik Hk , M µ i ∂ t Hi = -K i Ēi + a ik ∈F I S ik Ēk + a ik ∈F B S ik Ēk , (3)
where the mass matrices M γ i (γ stands for ε or µ), the stiffness matrix K i and the interface (or flux) matrix S ik write

(M γ i ) jl = γ i κ i ϕ i j • ϕ il , (K i ) jl = 1 2 κ i ϕ i j • ∇ × ϕ il + ϕ il • ∇ × ϕ i j , (S ik ) jl = 1 2 a ik ϕ i j • (ϕ kl × n ik ). ( 4
)
These matrices involve volume and surface integrals which are evaluated numerically over the master element, κ r , using the mapping Ψ κ i . For straight-sided simplices and parallelepipedic elements, all integrals can be computed using affine mappings. Since the Jacobian of the affine mapping is constant, the matrices in Eq. ( 4) can be precomputed and stored for κ r in advance of the main calculation once and for all. Straight-sided convex quadrilateral and hexahedral elements are expressed by isoparametric bilinear and trilinear mappings [START_REF] Knabner | Conditions for the invertibility of the isoparametric mapping for hexahedral finite elements[END_REF][START_REF] Knupp | On the invertibility of the isoparametric map[END_REF]. These mappings are not linear and the Jacobian is not constant over a quadrilateral and hexahedron anymore. Thus, all metric terms relating to Ψ κ i are evaluated, and need to be stored for each node in the elements. Finally, we also allow for curved elements which are usually present inside T h (e.g., around a boundary layer or a curved material interface) or in the vicinity of ∂Ω. These are introduced in order to avoid geometrical errors and to achieve optimal convergence for high-order polynomials.

In our implementation we use quadratic interpolation of curved boundaries, and define isoparametric mappings based on that approximation of the surface [START_REF] Barrett | Multilinear jacobians for isoparametric planar elements[END_REF][START_REF] Branets | Extension of a mesh quality metric for elements with a curved boundary edge or surface[END_REF][START_REF] Fahs | Improving accuracy of high-order discontinuous Galerkin method for time-domain electromagnetics on curvilinear domains[END_REF]. For such mappings, the Jacobians are not constants and all matrices should be computed and stored for each curved element.

Once the local mapping Ψ κ i is introduced for each κ i , the integrals in Eq. ( 4) are evaluated using cubature formulas which are able to accurately integrate high-order polynomials over the master element κ r . A cubature is a set of N c 2D or 3D points {λ c k } k=N c k=1 with N c associated weights {ω c k } k=N c k=1 , where the number of points N c depends on the desired maximum order of polynomial required to be accurately integrated. Surveys of cubature formulas were performed in [START_REF] Cools | Monomial cubature rules since "Stroud": a compilation-part 2[END_REF][START_REF] Cools | Monomial cubature rules since "Stroud": a compilation[END_REF] and an on-line database containing many of the best known rules is described in [START_REF] Cools | An encyclopaedia of cubature formulas[END_REF]. Here, we adopt the efficient high-order symmetric rules proposed in [16-18, 54, 56] for hypercubes and simplices. Finally, edge integrals are calculated using Gauss-Legendre quadrature. We evaluate

(M γ i ) jl = γ i k=N c k=1 ω c k L (ℓ i ) j L (ℓ i ) l J κ i λ c k , (K i ) jl = 1 2 k=N c k=1 ω c k L (ℓ i ) j (j κ i ) jl ∇ × L (ℓ i ) l + L (ℓ i ) l (j κ i ) jl ∇ × L (ℓ i ) j λ c k , (S ik ) jl = 1 2 k=N c k=1 ω c k L (ℓ i ) j L (ℓ i ) l ∇J κ i λ c k , (5)
where each term in the brackets of the right hand sides are evaluated at the cubature points, J κ i is the determinant of the Jacobian matrix, and j κ i = J κ i J -1 κ i . The set of local system of ordinary differential equations for each κ i , Eq. (3), can be formally transformed in a global system. To this end, we suppose that all electric (resp. magnetic) unknowns are gathered in a column vector E (resp. H) of size N g = n e i=1 N i . Then system (3) can be rewritten as

M ε ∂ t E = KH -AH -BH, M µ ∂ t H = -KE + AE -BE, (6)
where we have the following definitions and properties: (a) M ε , M µ and K are N g × N g block diagonal matrices with diagonal blocks equal to M ε i , M µ i , and K i respectively. M ε and M µ are symmetric positive definite matrices, and K is a symmetric matrix; (b) A is also a N g × N g symmetric block sparse matrix, whose non-zero blocks are equal to S ik when a ik ∈ F I ; (c) B is a N g × N g skew-symmetric block diagonal matrix, whose non-zero blocks are equal to S ik when a ik ∈ F B . Let S = K -A -B; the system (6) rewrites as

M ε ∂ t E = SH, M µ ∂ t H = -S † E. ( 7
)
where the symbol † stands for the matrix transpose. To define a fully discrete scheme, we divide the time interval [0, T ] into M t uniform subintervals such that 0 = t 0 < t 1 < • • • < t M t = T , where t n = n∆t and ∆t is the fixed time step. Time integration of the semi-discrete system, Eq. ( 7), is done using a fourth-order explicit leap-frog scheme. Let us denote E n = E(t n ) and H n+1/2 = H(t n+1/2 ), the fully discrete central DGTD method writes

               M ε E n+1 -E n ∆t = GH n+ 1 2 , M µ H n+ 3 2 -H n+ 1 2 ∆t = -G † E n+1 , (8)
where G = S I -∆t 2 24 (M µ ) -1 S † (M ε ) -1 S and I is the identity matrix. The resulting DGTD method is analyzed in [START_REF] Fahs | High-order leap-frog based discontinuous Galerkin method for the time-domain Maxwell equations on non-conforming simplicial meshes[END_REF][START_REF] Fahs | A high-order non-conforming discontinuous Galerkin method for time-domain electromagnetics[END_REF] on purely simplicial meshes where it is shown that the method is non-dissipative, conserves a discrete form of the EM energy and is stable under the CFL-like condition

∆t ≤ 2/ (M ε ) -1 2 G(M µ ) -1 2 ,
where • denotes any matrix norm, and the matrix (M γ ) -1 2 is the inverse square root of M γ (γ stands for ε or µ).

3. DGTD method on hybrid tetrahedral-hexahedral meshes

The numerical scheme

The set of elements κ i of the mesh is now assumed to be partitioned into two subsets: one made of the tetrahedral (triangular in 2D) elements and the other one gathering the hexahedral (quadrilateral in 2D) elements. In the following, these two subsets are respectively referred as S τ and S h . There is no need of a particular assumption on the connectivity of the two subsets. In order to further describe this scheme, we introduce additional definitions. First, the problem unknowns are reordered as

E = E τ E h and H = H τ H h ,
where subvectors with a τ subscript (resp., a h subscript) are associated to the elements of the set S τ (resp., the set S h ). We deduce from this partitioning of the unknown vectors the following decompositions of the system matrices

M ε = M ε τ 0 0 M ε h , M µ = M µ τ 0 0 M µ h , K = K τ 0 0 K h , B = B τ 0 0 B h .
where M ε τ/h and M µ τ/h are symmetric positive definite matrices, K τ/h are symmetric matrices and B τ/h are skewsymmetric matrices. The matrix A which involves the interface matrices S ik is decomposed as

A = A ττ A τh A hτ A hh ,
where A ττ and A hh are symmetric matrices, and A τh = A † hτ . Introducing the two matrices S τ = K τ -A ττ -B τ and S h = K h -A hh -B h , the global system of ordinary differential equations ( 7) can be split into two systems

       M ε τ ∂ t E τ = S τ H τ -A τh H h , M µ τ ∂ t H τ = -S † τ E τ + A τh E h ,        M ε h ∂ t E h = S h H h -A hτ H τ , M µ h ∂ t H h = -S † h E h + A hτ E τ . (9)
Then, the proposed DGTD scheme on hybrid meshes can be written as

                   M ε τ E n+1 τ -E n τ ∆t = G τ H n+ 1 2 τ -G τh H n+ 1 2 h , M µ τ         H n+ 3 2 τ -H n+ 1 2 τ ∆t         = -G † τ E n+1 τ + G τh E n+1 h , (10a)
                     M ε h       E n+1 h -E n h ∆t       = G h H n+ 1 2 h -G hτ H n+ 1 2 τ , M µ h          H n+ 3 2 h -H n+ 1 2 h ∆t          = -G † h E n+1 h + G hτ E n+1 τ , (10b)
where

G τ = S τ I - ∆t 2 24 (M µ τ ) -1 S † τ (M ε τ ) -1 S τ , G τh = A τh I - ∆t 2 24 (M µ h ) -1 A † τh (M ε τ ) -1 A τh , G h = S h I - ∆t 2 24 (M µ h ) -1 S † h (M ε h ) -1 S h , G hτ = G † τh .
The DGTD method on triangular-quadrilateral meshes or tetrahedral-hexahedral meshes, Eq. [START_REF]Discontinuous Galerkin Methods[END_REF], is referred to as DGTD-P p Q q in the following. If the mesh contains fully triangular or tetrahedral elements (resp. quadrilateral or hexahedral elements), then Eq. (10a) (resp. Eq. (10b)) is referred to as DGTD-P p method (resp. DGTD-Q q method).

Energy conservation

We define the following discrete electromagnetic energy

E n = 1 2 (E n τ ) † M ε τ E n τ + (H n-1 2 τ ) † M µ τ H n+ 1 2 τ + (E n h ) † M ε h E n h + (H n-1 2 h ) † M µ h H n+ 1 2 h . (11)
In the following, we shall verify that E n is exactly conserved in absence of absorbing boundary conditions, i.e., ∆E := E n+1 -E n = 0.

Lemma 1. Using the hybrid DGTD-P p Q q method, Eq. ( 10), for solving the Maxwell system, Eq. ( 1), with PEC boundary conditions only, the energy E n is exactly conserved.

Proof. We denote by E

n+ 1 2 τ/h = (E n τ/h + E n+1 τ/h )/2.
The variation of the energy E n through one time step is given by

∆E = (E n+ 1 2 τ ) † M ε τ [E n+1 τ -E n τ ] + 1 2 (H n+ 1 2 τ ) † M µ τ [H n+ 3 2 τ -H n-1 2 τ ] + (E n+ 1 2 h ) † M ε h [E n+1 h -E n h ] + 1 2 (H n+ 1 2 h ) † M µ h [H n+ 3 2 h -H n-1 2 h ] = ∆t(E n+ 1 2 τ ) † G τ H n+ 1 2 τ -∆t(E n+ 1 2 τ ) † G τh H n+ 1 2 h -∆t(H n+ 1 2 τ ) † G † τ E n+ 1 2 τ + ∆t(H n+ 1 2 τ ) † G τh E n+ 1 2 h + ∆t(E n+ 1 2 h ) † G h H n+ 1 2 h -∆t(E n+ 1 2 h ) † G hτ H n+ 1 2 τ -∆t(H n+ 1 2 h ) † G † h E n+ 1 2 h + ∆t(H n+ 1 2 h ) † G hτ E n+ 1 2 τ = ∆t(E n+ 1 2 τ ) † [G τ -G τ ]H n+ 1 2 τ -∆t(E n+ 1 2 τ ) † [G τh -G † hτ ]H n+ 1 2 h + ∆t(E n+ 1 2 h ) † [G h -G h ]H n+ 1 2 h -∆t(E n+ 1 2 h ) † [G hτ -G † τh ]H n+ 1 2 τ ,
which completes the proof.

Stability analysis

We shall prove that the electromagnetic energy is a positive definite quadratic form of all unknowns under a CFL-like condition on the time step ∆t.

Lemma 2. For the hybrid DGTD-P p Q k method, Eq. ( 10), the electromagnetic energy, E n , is a positive definite quadratic form of all numerical unknowns provided that

∆t ≤ 2 min 1 β τ + β τh , 1 β h + β τh with              β τ = (M ε τ ) -1 2 G τ (M µ τ ) -1 2 , β h = (M ε h ) -1 2 G h (M µ h ) -1 2 , β τh = (M µ τ ) -1 2 G τh (M ε h ) -1 2 , ( 12
)
where . denotes any matrix norm, and the matrix (M

γ τ/h ) -1 2 is the inverse square root of M γ τ/h (γ stands for ε or µ). Proof.
Using the second relations of Eq. (10a) and Eq. (10b), the energy E n , Eq. ( 11) can be written as

E n = 1 2 (E n τ ) † M ε τ E n τ + (H n-1 2 τ ) † M µ τ H n-1 2 τ -∆t(H n-1 2 τ ) † G † τ E n τ -∆t(H n-1 2 τ ) † G τh E n h +(E n h ) † M ε h E n h + (H n-1 2 h ) † M µ h H n-1 2 h -∆t(H n-1 2 h ) † G † h E n h -∆t(H n-1 2 h ) † G hτ E n τ .
The mass matrices M ε τ/h and M µ τ/h are symmetric positive definite and we can construct their square root (also symmetric positive definite) denoted by (M ε τ/h ) 1 2 and (M µ τ/h ) 1 2 , respectively. We deduce from Eq. ( 11) that

E n ≥ 1 2 (M ε τ ) 1 2 E n τ 2 + (M µ τ ) 1 2 H n-1 2 τ 2 -∆tβ τ (M µ τ ) 1 2 H n-1 2 τ (M ε τ ) 1 2 E n τ -∆tβ τh (M µ τ ) 1 2 H n-1 2 τ (M ε h ) 1 2 E n h + (M ε h ) 1 2 E n h 2 + (M µ h ) 1 2 H n-1 2 h 2 -∆tβ h (M µ h ) 1 2 H n-1 2 h (M ε h ) 1 2 E n h -∆tβ hτ (M µ h ) 1 2 H n-1 2 h (M ε τ ) 1 2 E n τ ,
where β τ , β h and β τh are defined in Eq. ( 12), and

β hτ = (M µ h ) -1 2 G hτ (M ε τ ) -1 2 .
Then estimate ab ≤ (a 2 + b 2 )/2 leads to

E n ≥ 1 2 (M ε τ ) 1 2 E n τ 2 + (M µ τ ) 1 2 H n-1 2 τ 2 -β τ ∆t 2 (M µ τ ) 1 2 H n-1 2 τ 2 + (M ε τ ) 1 2 E n τ 2 -β τh ∆t 2 (M µ τ ) 1 2 H n-1 2 τ 2 + (M ε h ) 1 2 E n h 2 + (M ε h ) 1 2 E n h 2 + (M µ h ) 1 2 H n-1 2 h 2 -β h ∆t 2 (M µ h ) 1 2 H n-1 2 h 2 + (M ε h ) 1 2 E n h 2 -β hτ ∆t 2 (M µ h ) 1 2 H n-1 2 h 2 + (M ε τ ) 1 2 E n τ 2 .
We then sum up the lower bounds for the E n to obtain

E n ≥ 1 2 1 -∆t β τ + β hτ 2 (M ε τ ) 1 2 E n τ 2 + 1 2 1 -∆t β τ + β τh 2 (M µ τ ) 1 2 H n-1 2 τ 2 + 1 2 1 -∆t β h + β τh 2 (M ε h ) 1 2 E n h 2 + 1 2 1 -∆t β h + β hτ 2 (M µ h ) 1 2 H n-1 2 h 2 .
Since the parameters ε and µ are piecewise constant and G τh = G † hτ , we have that β hτ = β τh . Then, under the condition proposed in Lemma 2, the energy E n is a positive definite quadratic form of the numerical unknowns E n τ/h and H n-1 2 τ/h , and Eq. ( 12) states a sufficient condition for the stability of the hybrid DGTD-P p Q q method.

Implementation issues

For generality, we shall limit the discussions to the 3D case and regard the 2D case as a natural simplification.

Nonconforming hybrid mesh generation

A simple way for generating a non-conforming hybrid tetrahedral-hexahedral mesh without pyramidal or prismatic elements is based on element-type conversion. This indirect approach starts with the generation of a conforming or non-conforming mesh with only one type of elements (pure tetrahedra or pure hexahedra). It is partially converted to hexahedral or tetrahedral mesh by applying a conventional subdivision scheme. For instance, a tetrahedron can be split into four hexahedra [START_REF] Carey | Hexing the Tet[END_REF] while a hexahedron can be decomposed into five or six tetrahedra [START_REF] Apel | Transformation of hexaedral finite element meshes into tetrahedral meshes according to quality criteria[END_REF]. The resulting hybrid mesh is non-conformal where the hybrid interface between a tetrahedron and a hexahedron is either a triangle or a quadrilateral, see Figure 1. Although this approach avoids some difficulties, it rapidly increases the number of elements and tends to introduce badly shaped elements which may deteriorate their numerical properties. We will not pursue this approach further here, and will propose an alternative method.

Hybrid interface Hybrid interface

Figure 1: Examples of hybrid interfaces when the hybrid tetrahedral-hexahedral mesh is generated using the element-type conversion approach. As shown in the left, if a tetrahedron is converted to hexahedra, the hybrid interface is an entire quadrilateral face of the hexahedron, where three of its four vertices are hanging nodes (filled circles). When a hexahedron is converted to tetrahedra, the hybrid interface is an entire triangular face of the tetrahedron. In this case, none of the three vertices of the hybrid interface is a hanging node.

In this work, we consider a more general strategy which consists of the following steps: first, the computational domain Ω is decomposed into subdomains Ω i such that Ω = i Ω i . Then, an individual mesh for each subdomain Ω i is generated using a suitable meshing algorithm and the most appropriate element type. The resulting meshes should not share nodes, edges, faces, or cells. Finally, we combine these submeshes to obtain the non-conformal hybrid mesh of the whole domain Ω. This strategy of non-conforming hybrid mesh generation has several advantages over other approaches based on hybrid conforming mesh. For instance, it allows one (i) to avoid the intermediate layer of pyramids as transition elements between triangular and quadrilateral faces; (ii) to relax the node-to-node matching requirement at the interface between the meshes; (iii) to change or replace certain portions of the domain without changing the entire mesh. These features are particularly useful for parametric modelling of complex geometries where the object is built from several parts designed separately. Figure 2(a)-(b) illustrates examples of hybrid meshes generated by this strategy.

After combining the meshes, we need to find the hybrid non-conforming interfaces between triangular and quadrilateral faces in order to compute the flux matrices, see Figure 2(c). First, we locate the boundary faces of each subdomain mesh and store them in a list including their element and face location. Next, we use the fact that only the boundary faces of two subdomains that shared boundaries may participate in a partial face-face connection. Let us assume that two subdomains Ω i and Ω j share a boundary. For each of the boundary faces of the mesh of Ω i , we perform two tests: the first is an intersection test with every member of the boundary face list of the mesh of Ω j . This just indicates whether or not an intersection exists. The second is a find-intersection test that involves computing the set of intersections between triangular and quadrilateral faces. To perform these intersection tests, we use the class Intersection provided by the Eberly's Wild Magic library [START_REF] Eberly | Game Physics[END_REF]. For the reader's convenience, we briefly describe the general idea of this class. The test-intersection is based on the method of separating axes [START_REF] Gottschalk | OBB-tree: a hierarchical structure for rapid interference detection[END_REF], a method for determining whether or not two convex objects are intersecting. A test for nonintersection of two convex objects is simply stated: if there exists a line for which the intervals of projection of the two objects onto that line do not intersect, then the objects do not intersect. Such a line is called a separating axis. For a pair of convex polygons in 3D, only a finite set of direction vectors needs to be considered for separation tests. The polygons normals are potential separating directions. If either normal direction separates the polygons, then no intersection occurs. However, if neither normal direction separates the polygons, then two possibilities exist for the remaining potential separating directions. The first case is that the polygons are coplanar and the remaining potential separating directions are those vectors in the common plane and perpendicular to the polygon edges. The second case is that the polygon planes are not parallel and that each polygon is split by the plane of the other polygon. The remaining potential separating directions are the cross products of two edges, one from each polygon. In summary, for the intersection between a triangular and a quadrilateral face in 3D, a maximum of 19 separating directions need to be tested in the worst case. Given a line, it is straightforward to project the vertices of a polygon onto the line and compute a bounding interval for those projections. The two intervals obtained from the two polygons are easily compared for overlap. If an intersection exist, we perform a find-intersection query using Boolean operations on polygons [START_REF] Hubbard | Constructive solid geometry for triangulated polyhedra[END_REF][START_REF] Rivero | Boolean operations on general planar polygons[END_REF]. Here the intersection between a triangular and a quadrilateral face (both lies in the same plane) is a convex polygon with at most seven vertices. Let Q be a quadrilateral and T be a triangle with vertices (v 1 , v 2 , v 3 ) that are counterclockwise ordered. We define half-planes, π i j , 1 ≤ i < j ≤ 3, bounded by the lines (v i v j ) and containing the third vertex of the triangle T . The intersection between T and Q is defined as

T ∩ Q = ((Q ∩ π 12 ) ∩ π 23 ) ∩ π 13 .
The problem is reduced to find the intersection between a polygon and a half-plane which is relatively straightforward. For a more complete exposition of the theory of intersection tests, we refer the reader to the books [48, Chap. Figure 2: Examples of non-conformal hybrid meshes generated using the combining meshes strategy. In (a) we show a triangular-quadrilateral mesh where the filled circles stand for hanging nodes. A 2D hybrid interface is a line segment with at least one of its end-points is a hanging node.

In (b) we show a 3D tetrahedral-hexahedral mesh. An example of 3D hybrid interfaces is illustrated in (c) where we show the intersection between a tetrahedron and four hexahedra. Each hybrid interface is a convex polygon with at most seven vertices, where at least one of them is a hanging node (filled circle).

Non-conforming hybrid interface matrix calculation

A key ingredient for the implementation of the DG method on hybrid meshes is the computation of the matrix A τh involving integrals over a hybrid interface. We concentrate here on the evaluation of such integrals. In fact, the other matrices in Eq. ( 9) which involve integrals over an element or an interface between two elements of the same type are evaluated as specified in Eq. ( 5). We just have to utilize the correct set of basis functions corresponding to the particular element type.

Let a ik = τ i ∩ h k be a hybrid interface between a tetrahedron τ i and a hexahedron h k . Here a ik is a n-gon with n = 3, 4, 5, 6, 7, that is, a ik is a convex polygon with n edges (see Figure 3). Let also R be a regular polygon in R 2 and Υ : R ∋ ξ ξ ξ → x ∈ a ik ⊂ R 3 be an explicit parametrization of a ik . The matrix A τh is block diagonal, whose non-zero blocks are equal to the interface matrix S ik when a ik is a hybrid interface. The matrix S ik is defined as

(S ik ) jl = a ik θ i j (x)ϑ kl (x)dx = R L (p i ) j • Ψ -1 τ i • Υ(ξ ξ ξ) L (q k ) l • Ψ -1 h k • Υ(ξ ξ ξ) ∇Υ(ξ ξ ξ) dξ ξ ξ (13)
where θ i j ∈ P p i (τ i ) and ϑ kl ∈ Q q k (h k ) are the basis functions of degrees p i and q k in the elements τ i and h k , respectively. The integrals in Eq. ( 13) are calculated numerically and stored for each hybrid interface a ik . To this end, two approaches are possible 1. Partitioning the interface a ik (a n-gon) into n triangles and then performing cubature rules over each triangle [START_REF] Dunavant | High degree efficient symmetrical Gaussian quadrature rules for the triangle[END_REF][START_REF] Wandzura | Symmetric quadrature rules on a triangle[END_REF]. For n-gons with n ≥ 4, we use the centroid of a ik to partition it into triangles. The matrix S ik is computed with the following formula

(S ik ) jl = a ik θ i j (x)ϑ kl (x)dx = s=n s=1 ∆ s θ i j (x)ϑ kl (x)dx, (14)
where ∆ s is a triangle formed by two successive vertices of a ik and its centroid as illustrated in Figure 4(a). The parametrization of ∆ s can be defined as x = Υ(ξ ξ ξ) = 3 i=1 φ i (ξ ξ ξ)v i , where ξ ξ ξ = (ξ 1 , ξ 2 ) ∈ T 2 , φ i are the barycentric coordinates in T 2 and v i denotes the three vertices of ∆ s . To compute the last integral in Eq. ( 14), we lay down N c cubature points λ n ∈ T 2 and weights ω n . This is achieved by

∆ s θ i j (x)ϑ kl (x)dx = N c n=1 L (p i ) j (λ τ n )L (q k ) l (λ h n )ω n ∇Υ(λ n ) , (15)
where

λ τ n = Ψ -1 τ i • Υ(λ n ) ∈ T 3 and λ h n = Ψ -1 h k • Υ(λ n ) ∈ D 3 .
2. An alternative is to use cubature rules over polygons. However, the design and development of such rules has not yet reached a mature stage. To our knowledge this approach has not been used previously in the context of DG methods. Recently, Mousavi et al. [START_REF] Mousavi | Generalized gaussian quadrature rules on arbitrary polygons[END_REF] presented an algorithm to build symmetric cubature rules for arbitrary regular polygons. This algorithm combines elements of group theory and numerical optimization and was originally proposed by Wandzura and Xiao [START_REF] Wandzura | Symmetric quadrature rules on a triangle[END_REF] for the construction of cubature rules over triangles. Let P n be a reference regular n-gon with vertices p m = (cos 2πm/n, sin 2πm/n) for m = 1, 2, . . . , n in the ξ ξ ξ-coordinate system, see Figure 4(b). The parametrization Υ from P n to a ik is obtained via the isoparametric mapping x = Υ(ξ ξ ξ) = n i=1 φ i (ξ ξ ξ)v i , where v i are the vertices of the n-gon a ik and φ i are the Laplace shape functions at point p ∈ P n defined as [START_REF] Sukumar | Conforming polygonal finite elements[END_REF]

φ i (ξ ξ ξ) = α i (ξ ξ ξ) n j=1 α j (ξ ξ ξ) , α j (ξ ξ ξ) = (1 -|ξ ξ ξ| 2 ) sin 3 (π/n) cos(π/n) A i-1 A i , ξ ξ ξ = (ξ 1 , ξ 2 ) ∈ P n ,
where α i (ξ ξ ξ) is the Laplace weight function, A i is the area of the triangle [p, p i , p i+1 ] and A i-1 is the area of the triangle [p, p i-1 , p i ]. The Laplace shape functions are non-negative, interpolate nodal data and form a partition of unity. They are also linear on the polygon boundaries which ensures that the parametrization Υ(ξ ξ ξ) is an affine mapping. Once Υ(ξ ξ ξ) is given explicitly, the integral in Eq. ( 13) is evaluated by laying down N c cubature points on P n and proceed as in Eq. ( 15).

We compare the two approaches by evaluating the integral of x 6 y 6 z 6 over the regular pentagon, hexagon and heptagon. The relative errors between exact and approximate values are shown in Figure 5 as a function of the number of cubature points. We observe that for a relative error of 10 -8 the polygonal cubatures require 3 to 5 times fewer points than the partitioning approach. Hence, for the numerical results given in Section 5, we adopt the second approach since it is more practical and requires the least number of cubature points for a given degree of accuracy.

Numerical results

Two-dimensional examples

We consider here the 2D TM z -polarized Maxwell equations, that is, we solve Eq. ( 1) for the fields H = (H x , H y , 0) and E = (0, 0, E z ). Three examples with analytical solutions are studied in detail to illustrate the accuracy and performance of the DG method on hybrid triangular-quadrilateral meshes. The exact time-domain solutions of these problems can be found in [START_REF] Fahs | Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic wave propagation[END_REF]. . . : Convergence curves of two numerical integration approaches used to evaluate the integral of x 6 y 6 z 6 over the regular pentagon (left), hexagon (center) and heptagon (right).

• • • • • • • • • • • • • • • • • • a a a ik ik ik • • • • a ik • • • a ik
ξ 2 T 2 ξ 1 ■ v 1 ∆ s (a) . . . . . . . . . . . . v v v v v v

Eigenmode in a PEC cavity

In the first example, we propose to validate the implementation of the proposed DGTD-P p Q q method on nonconforming hybrid meshes by performing a numerical convergence test. To this end, we consider the propagation of an eigenmode which is a standing wave of frequency f = 212 MHz and wavelength of 1.4 m in a unitary vacuum-filled square cavity with PEC walls. The computations are performed on three successively refined hybrid meshes, named MS i for i = 1, 2, 3. Each mesh MS i consists of 5 × 2 2i-1 triangles and 3 × 2 2i-1 quadrilaterals, as shown in Figure 6. We present results with interpolation orders, p and q, up to degree 8 such that p = q. For each interpolation order, we use the maximum time step, ∆t, allowed by the stability condition of the hybrid method. In Figure 7 (left) we illustrate the h-convergence graphs as a function of the total number of degrees of freedom (#DOF). As expected, for smooth solutions the convergence rate behaves as O(∆t 4 + h p ) even for non-conforming hybrid meshes. Also shown in Figure 7 (right) is the p-convergence rate in linear-log scale, observing an exponential decay of the error as exp(-s p p), with decay constant s p = log[E(p -1, q -1)/E(p, q)] ≃ 2, where E(p, q) is the L 2 error when the interpolation orders p and q are used. By inspecting the results in Figure 7, we observe that the gain in accuracy is notable when refining the mesh or increasing the interpolation orders. Furthermore, to achieve a certain accuracy level, the DG methods require less #DOF for higher interpolation orders, as typically found with other DG methods based on fully simplicial meshes [START_REF] Fahs | High-order leap-frog based discontinuous Galerkin method for the time-domain Maxwell equations on non-conforming simplicial meshes[END_REF][START_REF] Fahs | A high-order non-conforming discontinuous Galerkin method for time-domain electromagnetics[END_REF]. , for different values of the interpolation orders p and q with p = q. As expected, the h-convergence rate behaves as O(∆t 4 + h p ) for all orders. On the right, we plot the slope of the error for the mesh MS 1 , confirming the spectral convergence of the hybrid DGTD-P p Q q method.

Plane-wave scattering by a dielectric circular cylinder

As an illustration of the capability of the hybrid DG method to handle materials, let us consider plane wave scattering by a penetrable circular cylinder with a radius of 0.6 m consisting of a dielectric nonmagnetic material with a relative permittivity ε r = 12.0. The frequency of the incident wave is f = 300 MHz and the wavelength is λ = 0.288 m. The surrounding medium is assumed to be vacuum. The computational domain Ω is bounded by a square of side length 3.2 m centered at (0, 0). A Silver-Müller absorbing condition is applied on the boundary of the square. We compare the performance of the DGTD-P p Q q method using a non-conforming hybrid mesh with the DGTD-P p method on a fully triangular mesh, see Figure 8. The hybrid mesh is generated by combining submeshes of three subdomains, Ω 1 = [-0.35, 0.35] 2 m 2 (inside the cylinder), Ω 2 = [-0.7, 0.7] 2 m 2 \Ω 1 (close to cylinder and completely enclosing it) and

Ω 3 = [-1.6, 1.6] 2 m 2 \(Ω 1 ∪ Ω 2 ) such that Ω = Ω 1 ∪ Ω 2 ∪ Ω 3 .
The regions Ω 1 and Ω 2 are filled with regular quadrilaterals while Ω 3 is filled with unstructured triangles. The global hybrid mesh has 1888 triangles, 1204 quadrilaterals, 5044 conforming interfaces, 320 non-conforming hybrid interfaces and 120 nonconforming interfaces between quadrilaterals. The average edge length of this mesh is of 0.16λ. The triangular mesh is obtained by replacing the quadrilateral elements of the hybrid mesh by unstructured triangles. The resulting triangular mesh is conformal and consists of 4376 triangles and 6626 interfaces with an average edge length of 0.18λ. For the hybrid DGTD-P p Q q method, the interpolation order, q, is a pair q = (q 1 , q 2 ), where q 1 and q 2 are the interpolation orders used in the smallest and biggest quadrilaterals. The final simulation time has been set to T = 33.33 ns which corresponds to 10 periods of the incident wave. The results are listed in Table 1 in terms of accuracy and computational cost to reach the final simulation time for different p and q. Inspecting these results, we observe that for a fixed p, the hybrid method is 1.4 to 1.9 times more accurate and requires 1.2 to 2.1 times less computing time compared to the triangular method. On Figure 9 we compare the x-wise 1D distributions of the real part of the discrete Fourier transform of E z and H y components. We observe from this figure that, the main error of DG solutions arises at the material interfaces where the solution loses smoothness and the results confirm the ability of the hybrid DG method to accurately model problems with discontinuous coefficients and solutions. Finally, contour plots of the E z component at time t = T for numerical simulations performed with the DG method on the hybrid and the triangular meshes, are illustrated in Figure 10 for p = 2 and q = (2, 2). We observe very good agreement between the solutions obtained by both DG methods and the exact solution even with relatively low grid resolutions of six points per λ. The hybrid method shows slightly better results than the triangular method. Table 1: Comparison between the DGTD-P p Q q method on a hybrid triangular-quadrilateral mesh and the DGTD-P p method on a fully triangular mesh. For the hybrid method, the interpolation order q is a pair (q 1 , q 2 ) where q 1 and q 2 are the interpolation orders used in the smallest and biggest quadrilaterals, respectively. We give the maximum ∆t allowed by the stability condition, the L 2 -errors at time t = 33.33 ns, the #DOF and the computing time (CPU) in seconds to achieve the final time period.

P 2 Q (2,
2) (center) and triangular DGTD-P 2 (right) methods. Zoomed-in region to show more details. In all the plots, we use 13 equally spaced contours from E z = -1.5166 to E z = 1.5299.

Standing wave in a wedge-shaped PEC resonator

As an example with geometric singularities, we consider the propagation of a standing wave of frequency f = 1293 MHz and wavelength λ = 0.232 m inside a domain bounded by the curves y = tan(7π/4)x, x 2 + y 2 = 1/4 and the x-axis. The non-convex nature of the outer boundary renders some components of the solution highly singular, making it a challenging problem. We compare approximate solutions obtained using the DGTD-P p method formulated on triangular meshes with those resulting from the DGTD-P p Q q method formulated on hybrid non-conforming meshes. The triangular mesh consists of 6280 elements and 9531 interfaces with an average edge length of 0.057λ. The hybrid mesh contains 1745 triangles, 1820 quadrilaterals, 6228 conforming interfaces and 280 non-conforming interfaces. The global hybrid mesh has an average edge length of 0.082λ. Moreover, we can see that the hybrid mesh is locally refined in the vicinity of the re-entrant corner located at the origin. The resulting meshes are depicted in Figure 11. The simulation time has been set to T = 38.5 ns which corresponds to 50 periods. Results are summarized in Table 2 in terms of accuracy and computational costs. Here, the hybrid DG method produces numerical solutions that are one order of magnitude more accurate than those delivered by the triangular DG method. Figure 12 shows the contour lines of the H x component at t = T for numerical solutions obtained using the DGTD-P 4 and DGTD-P 2 Q 4 method. One can clearly observe a good qualitative agreement between the exact solution and the hybrid DGTD-P 2 Q 4 solutions. The triangular DGTD-P 4 solution produces some oscillations near the re-entrant corner due to bad resolution. This is confirmed in Figure 13 where we plot the x-wise 1D distributions of the H x and H y components along the line y = 5x + 0.025 passing nearby the corner singularity. Inspecting Table 2 and Figure 13, we observe that, for the DGTD-P p method, the error decreases slowly when increasing the order from p = 2 to p = 4. In contrast, the hybrid method offers faster convergence when increasing the order in the quadrilateral elements from q = 2 to q = 4. More specifically, the global error depends on the global regularity of the exact solution which motivates the use of lower polynomial order together with mesh refinement near singular points in order to recover standard convergence rates. We show the non-conforming, locally refined, hybrid triangular-quadrilateral mesh (left) and the conforming fully triangular mesh (right), used in the simulation of the wedge-shaped resonator. The hybrid mesh is refined near the re-entrant corner using triangular elements.

Table 2: Comparison between the hybrid DGTD-P p Q q method and the triangular DGTD-P p method. We give the maximum ∆t allowed by the stability condition, the L 2 -errors at time t = 38.5 ns, the #DOF and the computing time (CPU) in seconds to achieve the final time period.

Three-dimensional examples

Eigenmode in a spherical cavity

As a first verification of the general 3D framework, let us consider the propagation of the (0, 1, 1)-mode inside a PEC spherical cavity of unit radius. The resonant frequency is 131 MHz and the wavelength is λ = 2.29 m. The exact time-domain solution of this problem can be found in [START_REF] Balanis | Advanced Engineering Electromagnetics[END_REF]. We compare the accuracy and efficiency of the DG method using hexahedral, tetrahedral and hybrid meshes. First, a fully hexahedral mesh consisting of 6413 hexahedra and 19602 quadrilateral faces is generated. Then, the fully tetrahedral mesh is obtained by considering the hexahedral mesh and splitting each hexahedron to six tetrahedra. The resulting tetrahedral mesh has 38478 tetrahedra and 77682 triangular faces. Finally, the hybrid mesh contains 13068 tetrahedra, 216 regular hexahedra, 26136 triangular faces, 540 quadrilateral faces and 2916 hybrid non-conforming interfaces. These meshes are illustrated in Figure 14. The average edge lengths of the hexahedral, tetrahedral and hybrid meshes are respectively 0.0387λ, 0.0463λ and 0.0676λ which correspond to 26, 22 and 15 points per wavelength. We report on results obtained by the hexahedral DGTD-Q q , tetrahedral DGTD-P p and hybrid DGTD-P p Q q methods with interpolation orders p and q ranging from p = 1 up to p = 5 with p = q. The simulation time is fixed to T = 76 ns which corresponds to 10 periods. Table 3 summarizes the values of time steps used in the simulation. In Figure 15 (left) we show the L 2 error as a function of the interpolation order. Inspecting these results we observe several things. First, the gain in accuracy is notable when increasing the interpolation order for all methods. Here, the exponential decay of the error as exp(-2p) is observed for all cases. Second, the tetrahedral method is 3.9 to 7.9 times more accurate than the hexahedral method when increasing the interpolation order from p = 1 to p = 5. The results reveal that the accuracy of the hybrid method is almost identical to that of the tetrahedral method. Figure 15 (center) shows the ratio between the #DOF required by the hexahedral, tetrahedral and hybrid methods to achieve the final simulation time. The results indicate that the hybrid method is able to reduce the #DOF for given p and q. For instance, with p = q = 5, the L 2 error is 7.1 × 10 -5 for the hybrid mesh, 4.7 × 10 -4 for the hexahedral method and 6.0 × 10 -5 for the tetrahedral mesh. With the hybrid mesh, the #DOF is reduced by a factor of 1.8 compared to hexahedral mesh and by factor of 2.8 compared to tetrahedral mesh. Note that, the higher #DOF required by the tetrahedral mesh stems from the fact that the number of elements in the tetrahedral mesh is 6 times higher than those in the hexahedral mesh. Finally, the comparison in terms of computing times is displayed in Figure 15 (right). It can be observed that the reduction in the #DOF is translated into a corresponding reduction in computational cost. Here, the hybrid method allows for a reduction in computing time by a factor ranging from 2.4 to 8.0 compared to hexahedral method and by a factor of 3.5 compared to tetrahedral method. The lower computing times of the hybrid method partially stems from its larger time steps employed in hybrid mesh compared to those used in tetrahedral and hexahedral meshes.

Exposure of human head tissues to a localized source radiation

As a final example, we consider a more challenging problem which is concerned with the propagation of an EM wave emitted by a localized source in a realistic geometrical model of head tissues. Starting from magnetic resonance

X Y Z X Y Z X Y Z
Figure 14: Cut-away view of the hexahedral mesh (left), the tetrahedral mesh (center) and the hybrid tetrahedral-hexahedral mesh (right), used in the simulation of the spherical cavity.

Table 3: The time steps allowed by the stability conditions of the DG method using hexahedral, tetrahedral and hybrid meshes for different interpolation orders. In the left we show the error at the final simulation time as a function of the interpolation order, obtained using DG methods on tetrahedral (TET), hexahedral (HEX) and hybrid (HYB) meshes. The ratio between the number of degrees of freedom and the ratio between the computing times are illustrated in the center and right panels, respectively. images of the Visible Human 2.0 project [START_REF] Ratiu | Medicine Meets Virtual Reality 11 -NextMed: Health Horizon[END_REF], head tissues are segmented and the interfaces of a selected number of tissues are triangulated. Different strategies can be used in order to obtain a smooth and accurate segmentation of head tissues and interface triangulations as well. The strategy adopted in this example consists in using a variant of Chew's algorithm [START_REF] Chew | Guaranteed-quality mesh generation for curved surfaces[END_REF], based on Delaunay triangulation restricted to the interface, which allows to control the size and aspect ratio of interface triangles. Example of triangulations of the skin, skull and brain are shown in Figure 16. Then, these triangulated surfaces are used as inputs for the generation of volume meshes of the tissues. Finally, the GHS3D tetrahedral mesh generator [START_REF] George | Automatic mesh generator with specified boundary[END_REF] is used to mesh the volume domains between the various interfaces. The exterior of the head must also be meshed, up to a certain distance from the skin. The computational domain is here artificially bounded by a surface, denoted by Γ a , on which the absorbing boundary condition is imposed. Overall, the constructed geometrical model considered here consists of four tissues (skin, skull, CSF -Cerebro Spinal Fluid and brain). The characteristics of the tissues are summarized in Table 4 where the values of the relative permittivity ε r , and the conductivity σ, correspond to a frequency F = 1800 MHz and have been obtained from a special purpose online data base. For all tissues, the relative permeability, µ r , is set to 1. Given this tetrahedral mesh of the head tissues, the hybrid tetrahedral-hexahedral mesh is built in the following way (see the left panel of Figure 17): (a) the artificial boundary Γ a is defined as a plane-parallel surface; (b) a second plane-parallel surface, denoted by Γ 0 , lying close to the head and completely enclosing it, is introduced; (c) the volume space between Γ a and Γ 0 (vacuum space) is discretized using a regular hexahedral mesh; (d) the volume space between Γ 0 and the skin is discretized with unstructured tetrahedral elements. The global hybrid mesh consists of 197691 tetrahedra, 7000 hexahedra, 394397 triangular interfaces, 21900 quadrilateral interfaces and 5124 hybrid non-conforming interfaces. The minimum, maximum and average edge lengths of the hybrid mesh are equal to 1.12 mm, 43.14 mm and 17.26 mm, respectively. To compare the performance of the DG methods on hybrid and tetrahedral meshes, a fully tetrahedral mesh is generated by considering the hybrid mesh and replacing the hexahedral part by unstructured tetrahedral elements. We stress that the resulting tetrahedral mesh is conforming i.e., the same triangulated surface Γ 0 is used for the generation of tetrahedral elements inside and outside Γ 0 , as shown in Figure 17 (right). The global tetrahedral mesh consists of 350581 tetrahedra and 709395 interfaces. The minimum, maximum and average edge lengths of the tetrahedral mesh are equal to 0.49 mm, 44.32 mm and 13.98 mm, respectively. Finally, a dipolar type source is localized near the right ear of the head yielding a current of the form J z (x, t) = δ(xx d ) f (t), where f (t) is a sinusoidally varying temporal signal and x d is the localization point of the source. The physical simulation time has been fixed to T = 2.78 ns which corresponds to five periods of the temporal signal. A discrete Fourier transform of the components of the electric field is computed during the last period of the simulation. For this problem, we consider using the DGTD-P 2 method on the tetrahedral mesh and the DGTD-P 2 Q 2 and DGTD-P 2 Q 1 methods on the hybrid mesh. The time steps used in the simulations are ∆t = 2.728E-04 ns for the DGTD-P 2 method, ∆t = 6.655E-04 ns for the DGTD-P 2 Q 2 method and ∆t = 1.012E-03 ns for the DGTD-P 2 Q 1 method. These values correspond to the maximum ∆t allowed by the stability conditions of the DG methods. Contour lines of the modulus of the electric field, |E|, on selected planes of the skin for the approximate solutions resulting from the hybrid DGTD-P 2 Q 2 and tetrahedral DGTD-P 2 methods are visualized in Figure 18. Visually both methods provide almost the same solutions except slight amplitude differences around the right ear of the head. For the sake of completeness, we compare in Figure 19 the time evolution of the E z and H y components at a selected point in the center of the head. The comparisons show a visually perfect match between all solutions and no numerical artifacts due to the non-conformity of the mesh is observed. In Figure 19 we also plot the differences between the DGTD-P 2 solution and the hybrid solutions, showing that the errors are relatively very small (less than 0.1% for DGTD-P 2 Q 1 and less than 0.01% for DGTD-P 2 Q 2 ). Finally, we give some informations on the simulation times required by both schemes. On a workstation equipped with an Intel Xeon CPU E5-2650 @2.00 GHz and 64 GB of RAM, the DGTD-P 2 method requires 27 h 47 min for a total of 10183 time steps, the DGTD-P 2 Q 2 method requires 7 h 35 min for a total of 4177 time steps, and the DGTD-P 2 Q 1 method requires 4 h 17 min for a total of 2746 time steps. We obtain that in this case the hybrid DG method is 3.7 to 6.6 times less costly than the tetrahedral DG method.

p = q = 1 p = q = 2 p = q = 3 p = q = 4 p = q = 5

Concluding remarks

We have reported on our efforts to design a high-order accurate DG method for solving the time-domain Maxwell equations on non-conforming hybrid meshes. The proposed method combines a central numerical flux with a fourthorder leap-frog time integration. Within each mesh element, the EM field components are approximated using arbitrary high-order nodal polynomials. The geometrical flexibility is achieved by using unstructured triangular or tetrahedral meshes around complex geometric details and computational efficiency is then improved by using quadrilaterals or hexahedra to fill the remainder of the domain. The resulting DGTD method conserves a discrete form of the EM energy, and is stable under some CFL-like conditions. The key issues for the practical implementation of the proposed method are the generation of the hybrid meshes and the computation of flux matrices at the interfaces between hexahedral and tetrahedral faces. First, the mesh strategy adopted in this paper relies on the decomposition of the computational domain into non-overlapping subdomains, each one is meshed independently using the appropriate element type. Then, the hybrid interfaces are found using the method of separating axis and Boolean operations over polygons. This strategy demonstrates the benefits of relaxing the conformity requirement in the discretization process in view of facilitating the construction of meshes for complex propagation configurations, and improving the overall efficiency of the simulations. It is noteworthy that, most, if not all, of commercial or non-commercial meshgenerating softwares currently produce only conforming hybrid meshes by adding transition layers of pyramids or prisms. Second, the flux matrices are calculated using cubature rules over polygons. For a given accuracy, these rules are able to save more than 65% of the cubature points compared to other approaches. The proposed DGTD method is implemented in two and three space dimensions and tests on different wave propagation problems involving complex realistic geometries and heterogeneous media are presented. The spectral convergence for p-refinement method and the optimal convergence rate for h-refinement method have been verified. Results also show excellent agreement with exact solutions and no numerical artifacts caused by the non-conformity of the mesh are observed. Compared to DG methods on single element-type meshes, the DG method on hybrid meshes leads to clear reduction in computing time which is highly dependent on the interpolation order. In 2D, it is found that the hybrid mesh is able to produce at least a similar accuracy and allows the DG method to save up to 55% of the computing time and the number of degrees of freedom. In 3D, the gain becomes more significant, as the computational cost is reduced by up to 85%. A comparison between a fully hexahedral mesh and a fully tetrahedral mesh (obtained by converting each hexahedron to six tetrahedra) reveals that the tetrahedral method is always more accurate and more efficient than hexahedral method, excepting for linear interpolation orders. The parallelization of the proposed hybrid DG method and the comparison with DG methods on conforming hybrid meshes are left for future study. Figure 19: The top panels show the time evolution of H y (left) and E z (right) components at a selected point in the head. These solutions are obtained using the fully tetrahedral DGTD-P 2 method (solid line) and the hybrid DGTD-P 2 Q 2 (dashed line) and DGTD-P 2 Q 1 (dotted line) methods. In the bottom panels we plot the differences between the tetrahedral and hybrid DG solutions denoted as P 2 -P 2 Q 2 and P 2 -P 2 Q 1 .

  a h e d ro n c o n v e rt e d to 4 h e x a h e d ra

  7 & 11]-[21, Chap. 15].

Figure 3 :

 3 Figure 3: Five possible configurations for the hybrid interface a ik between a quadrilateral face and a triangular face in 3D.

Figure 4 :Figure 5

 45 Figure 4: Two numerical integration schemes used for the evaluation of integrals over a polygon a ik . In (a) we show a scheme based on the partition of a ik into triangles ∆ s . The scheme illustrated in (b) consists of using cubature rules over regular polygons.

Figure 6 :

 6 Figure 6: Sequences of three successively refined hybrid meshes used for assessing convergence of the eigenmode problem.

Figure 7 :

 7 Figure7: On the left is shown the L 2 -error for the EM field at time t = 10 ns as a function of the mesh resolution h = (DOF) 1/2 , for different values of the interpolation orders p and q with p = q. As expected, the h-convergence rate behaves as O(∆t 4 + h p ) for all orders. On the right, we plot the slope of the error for the mesh MS 1 , confirming the spectral convergence of the hybrid DGTD-P p Q q method.

Figure 8 :

 8 Figure 8: We show the non-conforming hybrid triangular-quadrilateral mesh (left) and the conforming fully triangular mesh (right), used for computing scattering by a dielectric cylinder of size 0.6 m.

Figure 9 :Figure 10 :

 910 Figure9: Plots of the E z (left) and H y (right) components along the axis y = 0 at time t = 33.33 ns for the scattering problem. We compare the exact solution (solid line) with solutions obtained using the hybrid DGTD-P 3 Q (3,2) (square) and triangular DGTD-P 3 (cross) methods. The vertical dotted lines at x = ±0.6 represent the material interfaces which distinguish the dielectric region, |x| ≤ 0.6, with ε r = 12.0 from the vacuum regions, 0.6 < |x| < 1.6, with ε r = 1. On the right figure we show zoomed-in views illustrating that main error arises at the material interfaces where the solution loses smoothness.

Figure 11 :

 11 Figure 11: We show the non-conforming, locally refined, hybrid triangular-quadrilateral mesh (left) and the conforming fully triangular mesh (right), used in the simulation of the wedge-shaped resonator. The hybrid mesh is refined near the re-entrant corner using triangular elements.

Figure 12 :Figure 13 :

 1213 Figure 12: Contour lines of the H x component at time t = 38.5 ns for the exact solution (left) and for solutions resulting from the hybrid DGTD-P 2 Q 4 (center) and triangular DGTD-P 4 (right) methods. Zoomed-in region to show more details. In all the plots, we use 20 equally spaced contours from H x = -0.440 to H x = 0.402.

Figure 15 :

 15 Figure 15: In the left we show the error at the final simulation time as a function of the interpolation order, obtained using DG methods on tetrahedral (TET), hexahedral (HEX) and hybrid (HYB) meshes. The ratio between the number of degrees of freedom and the ratio between the computing times are illustrated in the center and right panels, respectively.

Figure 16 :

 16 Figure 16: Surface mesh of the skin (left), the skull (center) and the brain (right).

Figure 17 :

 17 Figure 17: Planar views in a selected plane of the non-conforming hybrid tetrahedral-hexahedral mesh (left) and the conforming fully tetrahedral mesh (right) used for modelling the head tissues.

Figure 18 :

 18 Figure 18: Contour lines of the modulus of the electric field, |E|, in selected cut planes of the head, using the hybrid DGTD-P 2 Q 2 method (left) and the tetrahedral DGTD-P 2 method (right).

Table 4 :

 4 Electromagnetic characteristics of the selected head tissues at frequency f = 1800 MHz.

	Tissue	ε r	σ (S/m)	wavelength (mm)
	Skin	38.87	1.18	26.73
	Skull	15.56	0.43	42.25
	CSF	67.20	2.92	20.33
	Brain	43.55	1.15	25.26

Acknowledgements

The author would like to thank Dr. Seyed Mousavi of the University of California at Davis, USA, for providing the quadrature points over polygons to calculate flux matrices at hybrid interfaces. We also thank Dr. Nicolas Roquet for several discussions concerning the intersection tests for hybrid mesh generation. Finally, the author acknowledge partial support from the French National Agency for Research (ANR-10-JCJC-0905).


Samuel Cruz-Lara

Patrice Bonhomme
email: patrice.bonhomme@loria.fr

Christophe De Saint-Rat

Laurent Romary
email: laurent.romary@loria.fr

A GENERAL XML-BASED DISTRIBUTED SOFTWARE ARCHITECTURE FOR ACCESSING AND SHARING RESOURCES

This paper presents a general xml-based distributed software architecture with the aim of accessing and sharing resources in an opened client/server environment. The paper is organized as follows: First, we introduce the idea of a "General Distributed Software Architecture". Second, we describe the general framework in which this architecture is used. Third, we describe the process of information exchange and we introduce some technical issues involved in the implementation of the proposed architecture. Finally, we present some projects which are currently using the proposed architecture in the domain of textual resources.

Introduction.

) through which it might be possible to access and to share resources, which would be spread among different servers.

Figure 1. "General Architecture of the Network".

There are different reasons why to go about in this way. First, it should be considered that there already exist several sites all around which actually act as brokers for different types of resources which are stocked in some kind of databases. The general idea is that one should be able to bring together all the existing resources, in order to provide any user with a global access to these. In some cases this will lead to more coherence between these databases since redundancies might be detected or complementary resources put into correspondence. Second, the classical view of a centralised database containing all the information in a given domain is far from applicable to most kind of resources, because of their intrinsic diversity, for example. Also, there is a need for these resources to be created, and above all, maintained at a place where there is the competence to do so. For instance, textual resources may correspond to a wide variety of languages, genres and temporal periods which which are more easily maintained by groups of scholars possessing the proper competence. Finally, there can be specific constraints that can preclude some given resources to be deported to another site than the site which has originally created them. For example, we can have strong conditions related to property rights on the actual distribution of electronic versions of some textual resources. It is thus more sensible not to take the risk of hampering the agreement by overly spreading the corresponding contents. In the proposed network, each resource is only accessible through specific queries which can thus be controlled as to their actual applicability. Beyond that, we have also tried to base our proposals and/or implementations on existing (and even emerging) technologies (i.e. Java, CORBA, XML), so that to ensure some kind of durability to our work.

A General Framework

From the user's point of view, there should not be much change in the way the various resources funds are to be accessed, which means that whether there are one or several servers should be most transparent to him. As we will see, adopting a distributed framework, as opposed to the classical view of a centralised database, induces several specific problems for which this paper is trying to provide some plausible answers. Among those, we will have to deal specifically with the problem of broadcasting queries to different servers and conversely combining the corresponding result sets. As an example, statistics can only be dealt with in our distributed architecture if part of the computation is kept on the side of the remote servers and part is carried out locally (on the access server). All tasks dealing with the user interface should be concentrated on the client side, while the servers should accomplish searching and other computationally intensive operations. The network has the following characteristics: 1. each server is an autonomous unit containing its own data; 2. each server act as a "broker" and transmit, if needed, the request to other servers in the network which are known to it; 3. one given server is accessible to registered users through a general purpose Java-compatible web browser.

A General scenario

Any authorised user will be provided with an environment which will lead him along the following steps: 1. connection to a local server; 2. user identification; 3. choice of working servers. Given the list of available servers -accessible through the local server -in the network, together with their respective server profiles, the user will select those servers which may provide the proper resources or the proper services (tools) he wants to access or to use ; 4. selection of a subset of resources. Through an iterative process of requests to the selected servers, the user will build up a virtual subset (i.e. by way of pointers to individual resources) upon which he will actually work.

User's side

A user may interact with the network via a simple graphical user interface (in fact, only a general purpose Javacompatible web browser is needed). This interface is designed with a non-technical user in mind. However, technically advanced features should be available in an intuitive way. Although all system components will have their own interface, due to different functions, they must have the same "look and feel" (i.e. surface and behaviour). This reduces the time the user needs to become acquainted with the network and contributes to the aspect of simplicity. This user interface will be implemented as a client at the level of which little, not to say no, resource processing is to take place (notion of thin client). Basically, the interface will allow a user to make his different queries and will display result sets according to some specific style-sheets associated with these.

Connection

To be in line with the idea of a decentralised network, we have considered that a given user should only have to be registered at one given site and that no central user database should have to be set up.

Workspace

In order for a user to select the servers he wants to work with, he must interact with a "working space". The main purpose of the working space is to provide the user with a graphical interface that allows him to have access to the network and work on it.

The working space will offer the following functions:  Selection of the working servers;  Selection of working resources;  Manipulation of the selected resources;  Definition of the user's preferences.

Working Servers Selection

Before any query session (resource selection), a user has the possibility to select among a set of online servers, the server(s) he wants to work with. At any moment of the process, the user can edit his list of working servers and modify it.

Working Resources Selection

The principle is to restrict progressively the choice of resources in such a way that at the end of the selection process the user only keeps the subset he wants to work with. The user has the possibility to query either the whole network or his selected working servers. For that purpose, the user will have access to a friendly query interface for editing, modifying and sending his queries to his local server. Another more general way for the user to select his working resources is to browse through the whole set of referenced resources and to select the resources he is interested in.

Once a user has selected his working resources, he can keep them by using the "shopping basket" paradigm. As a matter of fact, one user can handle several shopping baskets, which he can save together with his workspace.

Manipulation of the selected resources.

As the main goal of this architecture is to give a user an access to a large set of resources, it is most important to provide a working space dedicated to the resources he wants to work with. The user may also use a set of tools allowing him to perform some kind of operations, queries for example, on the selected resources.

Server's side

As we have said, each server is an autonomous unit containing its own data. We have also considered that a given user should only have to be registered at one given site. From the point of view of the network, this results in making a distinction between:  "specialised servers", to manage users and resources, and  "broker servers", to share all existing resources. A given server can bear one or the other status, or even both depending on the actual implementation.

Users and Resources

Each server should manage not only its own resource database, but also its own users database. The users database contains both, the general user identification information, and all information related to workspaces.

Sharing Resources

Upon requests each server may act as a broker, that is, in order to share all the existing resources, it transmits requests to the other servers in the network. Each time a request is broadcast from an access server to a remote server, an identification tag (user id and authorisation level) should be transmitted in order to evaluate the applicability of the request.

A General Distributed Software Architecture

The network architecture is based on three major actors:  The Users (or Clients) which have been described in the preceding sections.  The Network Management Unit (NMU).

The NMU should be considered as the heart of the network: it allows to link all the servers which are connected to network. The NMU also maintains a database of information related to these servers (i.e. names, addresses, profiles, etc.).  The servers.

As mentioned before, this is a set of independent but associated specialised servers

The "Network Management Unit"

The Network Management Unit (NMU) is the only persistent link between the different servers affiliated to the network. It has thus the following functions (see Figure 2):  it maintains the list of servers affiliated to the network, together with a general profile for each of them (containing the languages dealt with by the server, the categories of document it may provide, etc). To this end, each server is described by its name, its address (URL), and a flag describing its current status;  it should be able to answer any query from each server belonging to the network when they have to know the list of available servers, the information attached to them and their current status;  it updates the local database of the specialised servers through the network, each time a change has been made to the network description or when something has changed on the part of a given server;  it is accessible through a simple form-based GUI, to allow the person responsible for the network (the "network master") to add a new server, modify the characteristics of a given server and/or disconnect a server from the network. At this point, no automatic processing will be made from the servers to NMU to ensure full control of the quality and coherence of the information ; Because of its central position in the architecture, the NMU is subject to several constraints, which are to condition its implementation. In particular, it has to be reliable as to its content since it is in charge of providing each server with the proper information to make it communicate with the others. There are three related but somehow different problems concerning clients, servers and the NMU: 1. we must allow servers to communicate with clients; 2. these servers have to communicate to each other; 3. we also have to implement a mechanism allowing these servers to communicate with the NMU. These points will be described in detail in the following sections.

Information exchange

In order to solve client/server and server/server communication, we have used Java's "servlets" and CORBA.

Servlets are the equivalent of CGI programs, but as they run in a single Java Virtual Machine (JVM), they don't have the per-request overhead of starting a new process per request. Servlets persist across multiple requests and hence can easily maintain states such as open database connections over these requests. In addition, there are all the usual advantages of Java : platform portability, safety (a buggy servlet won't bring down the server), enterprise integration (e.g. JDBC, CORBA, RMI), etc. Of course, servlets are associated to an HTTP server.

Servlets allow us to solve some security problems that we may have when using applets.

For instance, all information concerning users, will be stocked as XML documents. All these XML documents are created at the user's level by a Java applet. Because of security reasons, Java applets may not create files directly on the server. As a consequence, the XML document cannot be written on the server. Thus, the client's applet and the server communicate via a servlet. We should note that, the client's applet and the servlet allowing communication with the server, must be "physically" downloaded from the same URL. Obviously, the protocol used by the servlet in the framework of applet/server communication is HTTP. This communication is made by means of a request, in the very same way as a CGI script. That is, in order for the applet to retrieve the XML document concerning a user -this document contains the user's login and password -, we only have to make the request by using the servlet's related URL. Users database, which are in part administered via the user's interface, are composed by a set of XML documents. All these XML documents may be manipulated simply as files or as a real XML database.

Technical Issues

The XML level

All information flowing through the network will be encoded in XML. The encoding system will obviously be defined by DTDs. These DTDs allows to encode user's working spaces (WS), requests on data and meta-data (QL), results for these requests(RS), information related to users(UI) , and so on.

The MIME level

Because of efficiency, in particular in the framework of broadcasting requests from a server to another, not all needed information is presented in XML documents. In general on the Internet, data like sound, images, etc,. are encapsulated on a MIME layer. MIME allows to assign a data type (i.e., content-type), as well as to encapsulate heterogeneous data on a single data flow called MIME/Multipart: on a single channel, one single connection is enough in order to exchange multiple data.

The HTTP level

We use the HTTP protocol as the main communication protocol. Using HTTP is extremely simple (two single Java instructions allow to open an HTTP connection). However, http being a "stateless" protocol., each time a new page is loaded, the user is effectively disconnected from the server and it keeps no information allowing to know who was the user and what he was doing. The global scenario is: 1. Open the client/server connection; 2. send a request from the client to the server; 3. send a response from the server to the client; 4. close the client/server connection. Thus, even after logging into a site, each page accessed must pass the username and the password back to the server to verify the user's right to access the page. The client application (the browser) and the server application (the Web server) have no concept of local variables, local method calls, or objects. This problem is solved in part by using MIME/Multipart. On the server's side, all requests related to meta-data are analysed by means of a "local transaction". From a technical point of view, a Transaction is integrated within a Java's Servlet Session.

Query broadcast and result processing

Transaction

Every access to the resources database is made by means of a local transaction instance. In the framework of a transaction, all accesses to the resources database are considered to be persistent until the end of this Transaction.

A Transaction may be open or closed:

 open(), while receiving a request, for example, or during a database update operation;  close(), in general by a user's request or because a "timeout" operation.

We can also test any time, if the transaction is active by using the isopen() operation.

Managing an XML database is also done by means of a transaction. This transaction allows to control all current operations:  commit(), in order to execute an operation  abort(), to stop an operation.

Transaction/Query Servlet

Activating a Transaction by using a Java's Servlet allows to simulate the persistence of the client / server connection via a ServletSession provided by the JSDK 1 API. This persistence allows a user to progressively keep all results he is interested in. Thus, even if the evaluation of a request returns a high number of results, the user doesn't need to have a high storage or process capacity, in order to "navigate" inside the result's set.

Security aspects

Two main security threats can be identified in the proposed architecture:  Controlling access to the NMU and especially to the administration methods.  Securing the communications between servers. This basically means that we need to:  Identify trusted clients on the NMU;  Identify users on servers so only the authorised ones get access to protected resources;  Make sure that no user can connect as someone else, in order to get specific privileges;  Make sure that intercepting the data exchanged between servers is useless. Identification of users to the servers is implemented through a classic login and password procedure. Our purpose here is to present different methods that should be implemented and tested in order to secure the transactions taking place in the Network, on the one hand, between servers themselves, and on the other one hand, between servers and the NMU.

The methods presented here are all based on well known technologies largely used on the Internet, and available for systems implemented in the Java language.

Securing access to the NMU (Network Management Unit)

Some of the methods provided by the NMU IDL interface are used for administrating the NMU database, and so have to be protected in a way a fake client implemented by some "entrusted" people could not execute them, and thus trash the NMU database. Therefore, we have to implement a system allowing only trusted clients to call these methods. Two options are available for securing access to the NMU. The first one, the CORBA Security service, is implemented as a standard CORBA service. The second one is based on a password system.

Implementing the proposed architecture.

1 Java Servlet Development Kit, integrated within JDK 1.2

The SILFIDE Project[1].

SILFIDE is a tool for sharing, congenially and thoughtfully, a knowledge on different aspects of the French language. It consists in a network of data processing servers together with the necessary support.

The aim of SILFIDE is not to integrate the totality of the contents (corpora, glossaries, tools) of the available resources within an academic community, but to allow any researcher to be informed of the existence of such contents, to get a relatively precise idea of them and to be informed of the methods of access. In the case of resources which are widely used or which do not raise particular problems when accessed, SILFIDE will be able to propose the automatic transfer of the corresponding data. Thus, SILFIDE is a help provided to all laboratories of the French-speaking community and to those who are interested in the study or the automatic treatment of French language. In this respect, French has to be the main language of our project. On the one hand, most of the data available on SILFIDE will be in French or associated to equivalent data in French (in the case of a parallel corpus for instance). On the other hand, French will be the meta-language associated to the management of the resources either at the level of their documentation or at the level of the access interface to the corpus. However, a description of the server in other languages (English or German for instance) would be useful.

In the beginning, SILFIDE should be able to answer the following questions which may be raised by an user:  What are the available data ?  Where are they available and under which format ?  What are the conditions of access ?  How (and by whom) have these data been compiled ?  What is the validation degree of the Resources ?  What are the tools available to manipulate these resources ? Besides the function of access to linguistic resources, it might be interesting to propose directly accessible "on-line" tools to users who do not have access to an elaborate computer environment. Concordances can then be thought out for a set of selected texts together with elementary lexical statistics (frequencies, reduced deviations, etc.) Moreover, SILFIDE should be helpful by compiling (and possibly by documenting) the tools available in the field of textual resources manipulation. It can be a matter of encoding data but also libraries functions dedicated to normalised data. These different additional functions will have to be progressively integrated to the successive versions of the SILFIDE server.

The MLIS/ELAN Project[2].

Companies, research teams and individuals involved in language engineering or in undertakings such as translating, dictionary making or philology require LARGE corpora, lexica and similar electronic resources. On the other hand such linguistic databases have already been created for most European languages. It can be observed that the need exists, the product exists and yet the two have until now rarely been brought together.

This paradox can largely be explained by the fact that the language resources in question often cannot be accessed EASILY and usually cannot be exploited using STANDARDISED procedures.

Most existing European databases are not widely publicised, can only be accessed on site and are designed only for the immediate needs of the owner or producer. As a result, the software is often neither user-friendly nor portable, with documentation and support being non-existent or scant. In addition, the lack of a clear pricing means that each potential user must negotiate individually with the provider.

Standardised procedures postulate that similar data is made available in a similar format and that a common query system can operate on it. Appreciable progress was made during the 90's regarding the standardisation of formats through projects such as the TEI. Concerning the second requirement nothing has been achieved so far.

Two non-profit international associations -PAROLE and TELRI -have joined forces to undertake the ELAN project in order to link the existing resources with their potential users throughout Europe.

In order to serve the electronic multilingual resource market ELAN plans:

 to reinforce or, where necessary, create international standards by designing a common query language (ELAN-CQL) and by providing standardised resources for the following languages : Albanian, Belgian French, Belorussian, Bulgarian, Catalan, Croatian, Czech, Danish, Dutch, English, Estonian, Finnish, French, German, Greek, Hungarian, Irish, Italian, Latvian, Lithuanian, Polish, Portuguese, Romanian, Russian, Serbian, Slovakian, Slovene, Spanish, Swedish and Uzbek.

 to operate a user community network with active awareness-raising measures, a clear copyright policy, user support, e-mail user groups, etc.

One of the objectives of the MLIS/ELAN project is to define a proper software environment through which it might be possible to access and/or distribute linguistic resources, which would be spread among different servers.

There are different reasons why to go about in this fashion. First, it should be considered that there already exist several sites around the world (most of the European ones being represented in ELAN) which actually act as brokers for different types of linguistic resources. The idea emerging behind the sole period of the ELAN project is actually to be able to bring together the corresponding resources, in order to provide any user with a global access to these. In some cases this will lead to more coherence between these databases since redundancies might be detected or complementary resources (e.g. parallel texts) put into correspondence. Second, the classical view of a centralised database containing all the information in a given domain is far from applicable to linguistic resources where, because of their intrinsic diversity (prose, theatre, poetry, newspaper articles, dictionaries, historical documents etc.), there is a need for them to be created and above all maintained at a place where there is the competence to do so. In particular, the encoding background adopted within ELAN, that is SGML (and its subset XML), allows one to continuously enrich documents with specific linguistic annotations like part of speech (POS) tags or proper names (to quote only two possibilities).

Finally, there can be specific constraints that can preclude some given resources to be deported to another site than the site which has originally created them. In particular, some of the partners within ELAN have specific agreements with publishers, which express strong conditions on the actual distribution of electronic files. It is thus more sensible not to take the risk of hampering the agreement by overly spreading the corresponding contents. In the case of the ELAN network, each resource is only accessible through specific queries which can thus be controlled as to their actual applicability.

Information Exchange.

As indicated before, all information flowing throw the network is encoded in XML. In SILFIDE we have defined a "SILFIDE Interface Language" which is obviously specified by a DTD. This DTD is composed by a set of four modules defining working spaces (WS), requests on meta-data (QL), results for these requests (RS), and information related to users (UI). The MLIS/ELAN project use the same DTD but an extension mechanism adds a new DTD called "Elan Query" that specifies the ELAN "Common Query Language" Language.

The MLIS/ELAN network, based on the SILFIDE architecture, provides a simple but general way to interface any native database management system. We just have to specify an interface "ElanQueryHandler" for handling : 1. Connection and Transaction; 2. Query; 3. Result Set and Enumeration.

To "plug" the database into an ELAN server it is necessary to provide a driver that implements the "ElanQueryHandler" interface functions.

To optimise the transaction efficiency of the results and to avoid an overload of the network, we have implemented a "double cache" mechanism to store and send the result sets of the remote and local servers. The user has the possibility to specify the maximum size of the result set he can handle at a time. The size of the cache is set by the administrator of the server.

Figure 4 shows the information exchange process in the MLIS/ELAN network. We have the following steps : 1. A user send a CQL query to his local server; 2. The broadcast servlet intercepts and broadcast the query to each selected remote server (and even to the local server); 3. Each server establish a connection to the database, through the ElanQuery driver, and manage the transaction of the result with the server ("remote cache");

Figure 2 .

 2 Figure 2. "The Network Management Unit".

Figure 3 .

 3 Figure 3. "Information Flow in the Network".

<<

  !--SIL Document structure --> <!ELEMENT sil (uid, (%SIL.module.e;)+)> <!ATTLIST sil type (%SIL.x.a; workspace | user | query | resultset) "workspace" crdate CDATA #IMPLIED update CDATA #IMPLIED lang CDATA #IMPLIED sid CDATA #REQUIRED version CDATA #FIXED "0.5"> <!--lang : [RFC1766] language value sid : Silfide server identification --> <!--=================== --> <!--User IDentification --> <!--=================== --> <!--validated by the Silfide server admin. --> <!ELEMENT uid (login, passwd, access)> <!ATTLIST uid type (user | provider | both) "user"> default group IDREF #IMPLIED> <!--cross reference to a <group id="..."> element --> <!ELEMENT group EMPTY> <!ATTLIST group id ID #REQUIRED> <!--id : group identification key : inalf, lingua, etc. --> ]>