- Email: [email protected]
High-order discontinuous Galerkin method for time-domain electromagnetics on non-conforming hybrid meshes
Accepted Manuscript Title: High-order discontinuous Galerkin method for time-domain electromagnetics on non-conforming hybrid meshes Author: Hassan Fahs PII: DOI: Reference:
S0378-4754(14)00087-1 http://dx.doi.org/doi:10.1016/j.matcom.2014.04.002 MATCOM 4058
To appear in:
Mathematics and Computers in Simulation
Received date: Revised date: Accepted date:
19-8-2013 19-3-2014 7-4-2014
Please cite this article as: Hassan Fahs, High-order discontinuous Galerkin method for time-domain electromagnetics on non-conforming hybrid meshes, Mathematics and Computers in Simulation (2014), http://dx.doi.org/10.1016/j.matcom.2014.04.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript
High-order discontinuous Galerkin method for time-domain electromagnetics on non-conforming hybrid meshes Hassan Fahs∗
ip t
GPEM Group, Department of Materials & Structures, Nantes University, IFSTTAR, Route de Bouaye, 44344 Bouguenais Cedex, France
cr
Abstract
We present a high-order discontinuous Galerkin (DG) method for solving the time-dependent Maxwell equations
us
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
an
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
M
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
d
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
Ac ce pt e
methods allows for a notable reduction in computing time without sacrificing accuracy. Keywords: Maxwell’s equations, discontinuous Galerkin method, hybrid meshes, non-conformal meshes
1
1. Introduction
2
Nowadays, a variety of modeling strategies exist for the computer simulation of electromagnetic (EM) wave
3
propagations in the time-domain. In particular, the ability to deal with complex geometries has been extensively
4
studied in the context of unstructured methods as alternatives to the well established FDTD method [54]. In spite
5
of its versatility, the conventional FDTD method suffers from serious accuracy degradation when modeling curved
6
objects and small geometrical features. This is due to the commonly used staircase approximation (Cartesian grid)
7
which introduces large discretization errors, even when the grid size is very small. In contrast, unstructured solvers
8
such as the finite element time domain (FETD) [44] or the finite volume time domain (FVTD) [45] methods can
9
easily handle complex boundaries by dealing with unstructured tetrahedral meshes. Both methods are more expensive ∗ Corresponding
author Email address: [email protected], [email protected] (Hassan Fahs)
Preprint submitted to Elsevier
March 19, 2014
Page 1 of 28
than the FDTD method in terms of computing time and memory requirement. This motivates the development of
11
hybrid solution strategies, in which the computational efficiency of the Cartesian grid is combined with the geometric
12
flexibility of the unstructured mesh. The FDTD is then used in large homogeneous volumes of the computational
13
domain while the FVTD or FETD methods are employed in small regions near small scale features. For instance,
14
hybrid FETD-FDTD methods have been developed [56, 1], but they suffer from spurious reflections and late-time
15
instabilities, which limit their usefulness. A stable hybrid FETD-FDTD method has been described in [48]. To
16
achieve the time-stability, the scheme requires an implicit time integration in the unstructured elements but also a
17
transition layer of pyramids between the FDTD cubes and tetrahedra to ensure the mesh conformity. Low order
18
hybrid FDTD-FVTD methods are proposed in [23, 31]. These methods require the use of an intermediate layer of
19
hexahedra between the Cartesian grid and the tetrahedral mesh. Each cell in the layer is treated by the FVTD method,
20
which offers a correspondence between a quadrilateral face and two triangles in terms of the exchange of fluxes.
21
Again, weak instabilities have been reported for long-time simulations.
us
cr
ip t
10
The last few years have seen an increasing interest in so-called discontinuous Galerkin time domain (DGTD)
23
methods. The DG methods are a class of FE methods based on completely discontinuous piecewise polynomial
24
spaces for the numerical solution and the test functions. To obtain highly accurate and stable DG methods, suitable
25
numerical fluxes need to be designed over elemental interfaces. We refer to the lecture notes [10] and the textbook
26
[38] for details and history about DG methods. In particular, DG methods can easily handle elements of various
27
types and shapes [32], irregular non-conforming meshes [28, 15], and even locally varying polynomial degree [24],
28
and hence offer great flexibility in the mesh design. They also lead to (block-) diagonal mass matrices and therefore
29
yield fully explicit, inherently parallel methods when coupled with explicit time schemes. Furthermore, DG methods
30
are flexible with regards to the choice of the time-stepping scheme [27]. The DG methods have now acquired a
31
level of maturity and have successfully penetrated several scientific and technological communities following to their
32
adaptation to increasingly complex modeling contexts [29, 7, 50]. It is worthwhile to note that the DG method has
33
also been adopted for the first time in a commercial software as the time domain alternative of a very well known
34
EM wave simulation tool [52]. While DG methods on tetrahedral meshes have been extensively studied [37, 29, 15],
35
their development on hexahedral meshes has received relatively less attention [11]. One reason is that tetrahedral
36
meshing is highly automated while high quality hexahedral meshing of complex geometries commonly requires user
37
intervention and is labor intensive. This motivates the development of DG methods on conforming hybrid meshes
38
consisting of tetrahedra and hexahedra with transition layers of prisms or pyramids [5, 51]. Despite these advances,
39
the exploitation of non-conforming hybrid meshes in the context of a DG-based numerical methodology applied to
40
wave propagations, has been the subject of a very few studies so far. A 2D contribution in this context is the work
41
of Durochat et al. [19, 20] who proposed a DG method formulated on hybrid rectangular-triangular meshes for
42
Maxwell’s equations.
43
Ac ce pt e
d
M
an
22
We are concerned here with the design of a high-order nodal DGTD method for solving Maxwell’s equations 2
Page 2 of 28
on non-conforming hybrid meshes in 2D and 3D. The proposed method is based on centered numerical fluxes and
45
a fourth-order leap-frog time integration scheme. The hybrid mesh adopted here combines tetrahedral elements for
46
the discretization of irregularly shaped objects with hexahedral elements to fill the rest of the domain. Thus, the
47
resulting mesh is non-conformal in the sense that hanging nodes are allowed and no additional element types (e.g.,
48
pyramids or prisms) are needed to join tetrahedra with hexahedra. The non-conformal connection between these
49
elements is ensured thanks to the correct computation of flux terms by using efficient cubature rules over polygons.
50
The motivation behind this work stems partly from the following observations: first, the final mesh of an object is made
51
of separately constructed meshes and a conforming assembling is a very hard task. A good example to illustrate this
52
point is the interior penalty DGTD method formulated on non-conforming tetrahedral meshes presented in [15] for the
53
EM simulation of 3D integrated circuit packaging. For the complicated geometries involved in such applications, the
54
construction of a conforming mesh is a very hard task. Instead, a subdomain-based meshing approach is particularly
55
attractive and allowing meshes with different element-types is highly desirable in this context. Second, regarding the
56
computational effort, it turns out to be more efficient to compute on hexahedral meshes instead of tetrahedral meshes
57
to reach a desired error level. Therefore, the performance of the DG method can be increased by combining tetrahedra
58
with hexahedra. The outline of the paper is as follows. In Section 2 we briefly recall the basic features of DGTD
59
discretizations based on totally centered numerical fluxes and formulated on single element-type meshes composed
60
of triangles or quadrilaterals in 2D and tetrahedra or hexahedra in 3D. Then, we formulate the DG method on non-
61
conforming hybrid triangular-quadrilateral or tetrahedral-hexahedral meshes in Section 3. Subsequently, we analyze
62
the stability of the resulting algorithm and establish a sufficient CFL-like condition. Next, in Section 4, we discuss
63
some practical implementation issues including the generation of non-conforming hybrid meshes and the computation
64
of the flux matrices at hybrid interfaces between adjacent elements of different types. Finally, we present numerical
65
results for several test-cases in Section 5, and draw some conclusions of the study in Section 6.
66
2. Governing equations and DGTD formulation
68
cr
us
an
M
d
Ac ce pt e
67
ip t
44
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)
69
where ε and µ are respectively the permittivity and the permeability of the medium. These equations are posed on
70
a bounded domain Ω of Rd , d ≥ 2. At the domain boundary ∂Ω, we apply either perfect electric conductor (PEC)
71
conditions (n × E = 0) or absorbing boundary conditions (n × E + cµn × (n × H) = 0). Here n denotes the unit outward
72
normal to ∂Ω and c = (εµ)−1/2 is the speed of propagation.
73
We assume we dispose of a regular or irregular partition Th of Ω into ne non-overlapping elements κi , which can
74
be triangles or quadrilaterals in 2D and tetrahedra or hexahedra in 3D. Thus a (d −1)-dimensional face of each element 3
Page 3 of 28
75
76
77
78
κi in Th is allowed to contain hanging nodes. We shall suppose that each κi ∈ Th 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 ) P for all κi ∈ Th , where κr is either the open unit simplex T d = {ξξ = (ξ1 , . . . , ξd ) ∈ Rd : ξi ≥ 0, di=1 ξi ≤ 1} or the open unit hypercube Dd = {ξξ = (ξ1 , . . . , ξd ) ∈ Rd : 0 ≤ ξi ≤ 1} in Rd . 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 ) ∈ Nd0 , we define
80
spaces of polynomials of total degree ℓ ≥ 1 as follows
ip t
79
Qℓ (Dd ) = span{ξξα : 0 ≤ αi ≤ ℓ, 1 ≤ i ≤ d},
The dimension of each of these spaces depends on the spatial dimension d and on the degree ℓ and is given by d if κr = Dd , (ℓ + 1) N(ℓ, d) = (ℓ + d)!/ℓ!d! if κr = T d .
an
us
81
cr
α| ≤ ℓ}. Pℓ (T d ) = span{ξξα : 0 ≤ |α
Nℓ As a particular basis for Qℓ (Dd ) and Pℓ (T d ) consider the fundamental polynomials {L(ℓ) k }k=1 of multivariate Lagrange
83
interpolation, corresponding to a set of nodal points Sℓ = {η j , j = 1, . . . , Nℓ }, such that L(ℓ) k (η j ) = δ jk for η j ∈ Sℓ .
84
The distribution of the nodes η j is a key issue for the properties of the interpolation, especially for very high-order
85
approximations. For the unit simplex T d , we choose the nodal set derived from the electrostatics principle [35, 36],
86
while the Gauss-Legendre-Lobatto points are used for the hypercube Dd .
d
88
To each κi ∈ Th we assign an integer ℓi ≥ 1 that is the local polynomial degree inside κi ; collecting the ℓi and Ψκi in the vectors ℓ = {ℓi : κi ∈ Th } and Ψ = {Ψκi : κi ∈ Th }, respectively, we introduce the finite element space
Ac ce pt e
87
M
82
n d Vp (Ω, Th , Ψ ) = u ∈ [L2 (Ω)]d : u|κi ◦ Ψκi ∈ [Qℓi (κr )]d if Ψ−1 κi (κi ) = D ,
o d and u|κi ◦ Ψκi ∈ [Pℓi (κr )]d if Ψ−1 κi (κi ) = T , ∀κi ∈ Th ,
89
where [L2 (Ω)]d is the space of square-integrable functions on Ω. For two distinct elements κi and κk in Th , the
90
intersection aik = κi ∩ κk is a (d − 1)-dimensional face in Th , with unitary normal vector nik , oriented from κi towards
91
κk . For the boundary faces, the index k corresponds to a fictitious element outside the domain. Finally, we denote by
92
F I and F B the union of all interior and boundary faces of Th , respectively.
93
The derivation of the DG scheme for the simulation of electromagnetic wave propagation on purely simplicial
94
meshes has already been published in previous works [24, 30]. Therefore, we avoid the repetition of a comprehensive
95
derivation and refer the reader to these works for details. In brief, the electric and magnetic fields are approximated
96
inside each finite element κi by a linear combination of basis functions ϕi j = L(ℓj i ) ◦ Ψ−1 κi of degree ℓi with support κi
97
and with time-dependent coefficient functions Ei j (t) and Hi j (t) as follows Ei (x, t) =
X
Ei j (t)ϕi j (x),
Hi (x, t) =
X
Hi j (t)ϕi j (x).
(2)
1≤ j≤Ni
1≤ j≤Ni
4
Page 4 of 28
Here, the index j indicates the jth basis function and Ni = N(ℓi , d) is the number of degrees of freedom inside κi .
99
¯ i , respectively, the column vectors (Ei j )1≤ j≤Ni and (Hi j )1≤ j≤Ni . As usual for DG schemes, the We denote by E¯ i and H
100
Maxwell system, Eq. (1), is multiplied by a test function ϕ ∈ span{ϕi j , 1 ≤ j ≤ Ni } and integrated over each κi .
101
After integration by parts, inserting the DG approximation, Eq. (2), and after applying a centered numerical flux in
102
the boundary integrals, the semi-discrete formulation of the scheme in the physical element κi reads as =
¯i− Ki H
X
¯k− Sik H
X
¯k + Sik E
aik ∈F I
¯i Mµi ∂t H
= −Ki E¯ i +
¯ k, Sik H
X
Sik E¯ k ,
aik ∈F B
aik ∈F I
aik ∈F B
(3)
where the mass matrices Mγi (γ stands for ε or µ), the stiffness matrix Ki and the interface (or flux) matrix Sik write Z Z Z 1 1 ϕi j · ϕil , (Ki ) jl = (Mγi ) jl = γi ϕi j · ∇ × ϕil + ϕil · ∇ × ϕi j , (Sik ) jl = ϕi j · (ϕkl × nik ). (4) 2 κi 2 aik κi
us
103
X
cr
¯i Mεi ∂t E
ip t
98
These matrices involve volume and surface integrals which are evaluated numerically over the master element, κr ,
105
using the mapping Ψκi . For straight-sided simplices and parallelepipedic elements, all integrals can be computed using
106
affine mappings. Since the Jacobian of the affine mapping is constant, the matrices in Eq. (4) can be precomputed and
107
stored for κr in advance of the main calculation once and for all. Straight-sided convex quadrilateral and hexahedral
108
elements are expressed by isoparametric bilinear and trilinear mappings [41, 40]. These mappings are not linear and
109
the Jacobian is not constant over a quadrilateral and hexahedron anymore. Thus, all metric terms relating to Ψκi are
110
evaluated, and need to be stored for each node in the elements. Finally, we also allow for curved elements which are
111
usually present inside Th (e.g., around a boundary layer or a curved material interface) or in the vicinity of ∂Ω. These
112
are introduced in order to avoid geometrical errors and to achieve optimal convergence for high-order polynomials.
113
In our implementation we use quadratic interpolation of curved boundaries, and define isoparametric mappings based
114
on that approximation of the surface [26, 4, 6]. For such mappings, the Jacobians are not constants and all matrices
115
should be computed and stored for each curved element.
Ac ce pt e
d
M
an
104
116
Once the local mapping Ψκi is introduced for each κi , the integrals in Eq. (4) are evaluated using cubature formulas
117
which are able to accurately integrate high-order polynomials over the master element κr . A cubature is a set of N c
118
k=Nc k=Nc with N c associated weights {ωck }k=1 , where the number of points N c depends on the desired 2D or 3D points {λck }k=1
119
maximum order of polynomial required to be accurately integrated. Surveys of cubature formulas were performed
120
in [14, 12] and an on-line database containing many of the best known rules is described in [13]. Here, we adopt
121
the efficient high-order symmetric rules proposed in [42, 55, 17, 16, 18] for hypercubes and simplices. Finally, edge
5
Page 5 of 28
integrals are calculated using Gauss-Legendre quadrature. We evaluate γ
(Mi ) jl
=
γi
k=N Xc k=1
(Ki ) jl
=
(Sik ) jl
=
h i i) , ωck L(ℓj i ) L(ℓ kJ k κ i l c λk
k=N i 1 Xc c h (ℓi ) ωk L j (jκi ) jl ∇ × Ll(ℓi ) + Ll(ℓi ) (jκi ) jl ∇ × L(ℓj i ) c , λk 2 k=1 k=N i 1 Xc c h (ℓi ) (ℓi ) ωk L j Ll k∇Jκi k c , λk 2 k=1
(5)
ip t
122
where each term in the brackets of the right hand sides are evaluated at the cubature points, kJκi k is the determinant of
124
the Jacobian matrix, and jκi = kJκi kJ−1 κi .
127
us
126
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 Pe (resp. H) of size Ng = ni=1 Ni . Then system (3) can be rewritten as Mε ∂t E
=
Mµ ∂t H
= −KE + AE − BE,
KH − AH − BH,
an
125
cr
123
(6)
where we have the following definitions and properties: (a) Mε , Mµ and K are Ng × Ng block diagonal matrices with
129
diagonal blocks equal to Mεi , Mi , and Ki respectively. Mε and Mµ are symmetric positive definite matrices, and K is
130
a symmetric matrix; (b) A is also a Ng × Ng symmetric block sparse matrix, whose non-zero blocks are equal to Sik
131
when aik ∈ F I ; (c) B is a Ng × Ng skew-symmetric block diagonal matrix, whose non-zero blocks are equal to Sik
132
when aik ∈ F B . Let S = K − A − B; the system (6) rewrites as
Ac ce pt e
d
µ
M
128
Mε ∂t E
=
SH,
Mµ ∂t H
=
−S† E.
(7)
133
where the symbol † stands for the matrix transpose. To define a fully discrete scheme, we divide the time interval
134
[0, T ] into Mt uniform subintervals such that 0 = t0 < t1 < · · · < t Mt = T , where tn = n∆t and ∆t is the fixed time step.
135
Time integration of the semi-discrete system, Eq. (7), is done using a fourth-order explicit leap-frog scheme. Let us
136
137
denote En = E(tn ) and Hn+1/2 = H(tn+1/2 ), the fully discrete central DGTD method writes En+1 − En 1 Mε = GHn+ 2 , ∆t (8) 1 n+ 23 − Hn+ 2 † n+1 Mµ H = −G E , ∆t 2 µ −1 † ε −1 where G = S I − ∆t 24 (M ) S (M ) S and I is the identity matrix. The resulting DGTD method is analyzed in
138
[30, 25] on purely simplicial meshes where it is shown that the method is non-dissipative, conserves a discrete form
139
of the EM energy and is stable under the CFL-like condition 1
1
∆t ≤ 2/k(Mε )− 2 G(Mµ )− 2 k, 1
140
where k · k denotes any matrix norm, and the matrix (Mγ )− 2 is the inverse square root of Mγ (γ stands for ε or µ). 6
Page 6 of 28
141
3. DGTD method on hybrid tetrahedral-hexahedral meshes
142
3.1. The numerical scheme The set of elements κi of the mesh is now assumed to be partitioned into two subsets: one made of the tetrahe-
144
dral (triangular in 2D) elements and the other one gathering the hexahedral (quadrilateral in 2D) elements. In the
145
following, these two subsets are respectively referred as Sτ and Sh . There is no need of a particular assumption on
146
the connectivity of the two subsets. In order to further describe this scheme, we introduce additional definitions. First,
147
the problem unknowns are reordered as E τ E = E
us
h
where subvectors with a τ subscript (resp., a h subscript) are associated to the elements of the set Sτ (resp., the set Sh ). We deduce from this partitioning of the unknown vectors the following decompositions of the system matrices B K Mµ 0 Mε 0 0 0 τ τ τ τ . , B = , K = , Mµ = Mε = 0 Mε 0 B 0 K 0 Mµ h
h
151
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 Sik is decomposed as A ττ Aτh A = , A A
154
155
hh
where Aττ and Ahh are symmetric matrices, and Aτh = A†hτ . Introducing the two matrices Sτ = Kτ − Aττ − Bτ and
Ac ce pt e
153
d
hτ
152
h
h
M
150
an
149
cr
H τ and H = , H
h
148
ip t
143
Sh = Kh − Ahh − Bh , the global system of ordinary differential equations (7) can be split into two systems ε Mτ ∂t Eτ = Sτ Hτ − Aτh Hh , Mµτ ∂t Hτ = −S†τ Eτ + Aτh Eh , ε Mh ∂t Eh = Sh Hh − Ahτ Hτ , Mµh ∂t Hh = −S†h Eh + Ahτ Eτ . Then, the proposed DGTD scheme on hybrid meshes can be written as ! n+1 − Enτ n+ 1 n+ 1 ε Eτ M = Gτ Hτ 2 − Gτh Hh 2 , τ ∆t n+ 3 1 H 2 − Hn+ 2 µ τ τ = −G†τ En+1 + Gτh En+1 τ h , Mτ ∆t n+1 Eh − Enh ε = Mh ∆t n+ 3 1 H 2 − Hn+ 2 µ h h = M h ∆t
n+ 12
Gh Hh
n+ 21
− Ghτ Hτ
(9)
(10a)
, (10b)
−G†h En+1 h
+
Ghτ En+1 τ ,
7
Page 7 of 28
156
where
! ∆t2 µ −1 † ε −1 Gτ = Sτ I − (Mτ ) Sτ (Mτ ) Sτ , 24 ! ∆t2 µ −1 † ε −1 (Mh ) Sh (Mh ) Sh , Gh = Sh I − 24
! ∆t2 µ −1 † ε −1 I− (Mh ) Aτh (Mτ ) Aτh , 24
Gτh = Aτh Ghτ = G†τh .
The DGTD method on triangular-quadrilateral meshes or tetrahedral-hexahedral meshes, Eq. (10), is referred to as
158
DGTD-P p Qq in the following. If the mesh contains fully triangular or tetrahedral elements (resp. quadrilateral or
159
hexahedral elements), then Eq. (10a) (resp. Eq. (10b)) is referred to as DGTD-P p method (resp. DGTD-Qq method).
160
3.2. Energy conservation
cr
We define the following discrete electromagnetic energy En =
us
161
ip t
157
1 1 1 n † ε n n− 1 n− 1 µ n+ µ n+ (Eτ ) Mτ Eτ + (Hτ 2 )† Mτ Hτ 2 + (Enh )† Mεh Enh + (Hh 2 )† Mh Hh 2 . 2
(11)
In the following, we shall verify that En is exactly conserved in absence of absorbing boundary conditions, i.e.,
163
∆E := En+1 − En = 0.
164
Lemma 1. Using the hybrid DGTD-P p Qq method, Eq. (10), for solving the Maxwell system, Eq. (1), with PEC
165
boundary conditions only, the energy En is exactly conserved.
166
n Proof. We denote by Eτ/h2 = (Enτ/h + En+1 τ/h )/2. The variation of the energy E through one time step is given by
M
an
162
=
1 n+ 1 n+ 3 n− 1 µ ) Mετ [En+1 − Enτ ] + (Hτ 2 )† Mτ [Hτ 2 − Hτ 2 ] + τ 2 1 n+ 1 n+ 3 n− 1 n+ 21 † ε n+1 n (Eh ) Mh [Eh − Eh ] + (Hh 2 )† Mµh [Hh 2 − Hh 2 ] 2 n+ 12 †
(Eτ
Ac ce pt e
∆E
d
n+ 1
=
n+ 21 †
n+ 12
n+ 21
n+ 21
∆t(Eτ
) Gτ Hτ
∆t(Eh )† Gh Hh
=
n+ 21
n+ 1
n+ 12
) [Gτ − Gτ ]Hτ
∆t(Eh 2 )† [Gh − Gh ]Hh
167
which completes the proof.
168
3.3. Stability analysis
169
170
n+ 12
n+ 12
n+ 12
) Gτh Hh
− ∆t(Eh )† Ghτ Hτ
n+ 21 †
∆t(Eτ
n+ 21 †
− ∆t(Eτ
n+ 12 †
n+ 12
n+ 21
n+ 12
− ∆t(Hτ
) G†τ Eτ
− ∆t(Hh )† G†h Eh
n+ 12 †
n+ 21
n+ 1
n+ 1
− ∆t(Eτ
) [Gτh − G†hτ ]Hh
n+ 12 †
n+ 21
n+ 21
n+ 12
+ ∆t(Hτ
) Gτh Eh
+
+ ∆t(Hh )† Ghτ Eτ
+
− ∆t(Eh 2 )† [Ghτ − G†τh ]Hτ 2 ,
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.
8
Page 8 of 28
171
Lemma 2. For the hybrid DGTD-P p Qk method, Eq. (10), the electromagnetic energy, En , is a positive definite
172
quadratic form of all numerical unknowns provided that
1 1 , βτ + βτh βh + βτh
!
with
h
1
Proof. Using the second relations of Eq. (10a) and Eq. (10b), the energy En , Eq. (11) can be written as 1 n † ε n n− 1 n− 1 n− 1 n− 1 (Eτ ) Mτ Eτ + (Hτ 2 )† Mµτ Hτ 2 − ∆t(Hτ 2 )† G†τ Enτ − ∆t(Hτ 2 )† Gτh Enh En = 2 n− 1 n− 1 n− 1 n− 1 +(Enh )† Mεh Enh + (Hh 2 )† Mµh Hh 2 − ∆t(Hh 2 )† G†h Enh − ∆t(Hh 2 )† Ghτ Enτ .
us
174
where k.k denotes any matrix norm, and the matrix (Mγτ/h )− 2 is the inverse square root of Mγτ/h (γ stands for ε or µ).
cr
173
(12)
ip t
∆t ≤ 2 min
1 1 βτ = k(Mετ )− 2 Gτ (Mµτ )− 2 k, 1 µ 1 βh = k(Mεh )− 2 Gh (Mh )− 2 k, βτh = k(Mµτ )− 21 Gτh (Mε )− 21 k,
The mass matrices Mετ/h and Mµτ/h are symmetric positive definite and we can construct their square root (also sym-
176
metric positive definite) denoted by (Mετ/h ) 2 and (Mµτ/h ) 2 , respectively. We deduce from Eq. (11) that
1
≥
1
1 1 1 1 1 1 1 n− 1 n− 1 n− 1 k(Mετ ) 2 Enτ k2 + k(Mµτ ) 2 Hτ 2 k2 − ∆tβτ k(Mµτ ) 2 Hτ 2 k k(Mετ ) 2 Enτ k − ∆tβτh k(Mµτ ) 2 Hτ 2 k k(Mεh ) 2 Enh k 2 1 1 1 1 1 1 n− 1 n− 1 n− 1 +k(Mεh ) 2 Enh k2 + k(Mµh ) 2 Hh 2 k2 − ∆tβh k(Mµh ) 2 Hh 2 k k(Mεh ) 2 Enh k − ∆tβhτ k(Mµh ) 2 Hh 2 k k(Mετ ) 2 Enτ k ,
M
En
an
175
1
≥
d
En
178
1
where βτ , βh and βτh are defined in Eq. (12), and βhτ = k(Mµh )− 2 Ghτ (Mετ )− 2 k. Then estimate ab ≤ (a2 + b2 )/2 leads to 1 1 1 1 ∆t ∆t 1 1 1 µ 1 n− µ 1 n− µ 1 n− k(Mετ ) 2 Enτ k2 + k(Mτ ) 2 Hτ 2 k2 − βτ k(Mτ ) 2 Hτ 2 k2 + k(Mετ ) 2 Enτ k2 − βτh k(Mτ ) 2 Hτ 2 k2 + k(Mεh ) 2 Enh k2 2 2 2 1 1 ∆t ∆t 1 1 1 1 1 1 n− n− n− 1 +k(Mεh ) 2 Enh k2 + k(Mµh ) 2 Hh 2 k2 − βh k(Mµh ) 2 Hh 2 k2 + k(Mεh ) 2 Enh k2 − βhτ k(Mµh ) 2 Hh 2 k2 + k(Mετ ) 2 Enτ k2 . 2 2
Ac ce pt e
177
We then sum up the lower bounds for the En to obtain βτ + βhτ 1 1 1 − ∆t k(Mετ ) 2 Enτ k2 + En ≥ 2 2 βh + βτh 1 1 1 − ∆t k(Mεh ) 2 Enh k2 + 2 2
βτ + βτh 1 1 n− 1 1 − ∆t k(Mµτ ) 2 Hτ 2 k2 + 2 2 1 βh + βhτ 1 n− 1 1 − ∆t k(Mµh ) 2 Hh 2 k2 . 2 2
179
Since the parameters ε and µ are piecewise constant and Gτh = G†hτ , we have that βhτ = βτh . Then, under the condition
180
proposed in Lemma 2, the energy En is a positive definite quadratic form of the numerical unknowns Enτ/h and Hτ/h2 ,
181
and Eq. (12) states a sufficient condition for the stability of the hybrid DGTD-P p Qq method.
182
4. Implementation issues
183
n− 1
For generality, we shall limit the discussions to the 3D case and regard the 2D case as a natural simplification.
9
Page 9 of 28
184
4.1. Nonconforming hybrid mesh generation A simple way for generating a non-conforming hybrid tetrahedral-hexahedral mesh without pyramidal or prismatic
186
elements is based on element-type conversion. This indirect approach starts with the generation of a conforming or
187
non-conforming mesh with only one type of elements (pure tetrahedra or pure hexahedra). It is partially converted
188
to hexahedral or tetrahedral mesh by applying a conventional subdivision scheme. For instance, a tetrahedron can
189
be split into four hexahedra [8] while a hexahedron can be decomposed into five or six tetrahedra [2]. The resulting
190
hybrid mesh is non-conformal where the hybrid interface between a tetrahedron and a hexahedron is either a triangle
191
or a quadrilateral, see Figure 1. Although this approach avoids some difficulties, it rapidly increases the number of
192
elements and tends to introduce badly shaped elements which may deteriorate their numerical properties. We will not
193
pursue this approach further here, and will propose an alternative method.
us
M
Neighbor hexahedron
Hybrid interface
an
Hexahedron converted to 5 tetrahedra
Tet
ing
rm
.
fo on −c s on ace d n erf bri int
..
.
0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1111111111111 0000000000000 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 10 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 10 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111
Hy
rah e to 4dron c hex onv ahe erte dra d
Neighbor tetrahedron
Hybrid interface
111111 000000 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 000000 111111 000000 111111 0000000 1111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111
cr
ip t
185
Figure 1: Examples of hybrid interfaces when the hybrid tetrahedral-hexahedral mesh is generated using the element-type conversion approach. As
d
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
194
195
Ac ce pt e
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 S domain Ω is decomposed into subdomains Ωi such that Ω = i Ωi . Then, an individual mesh for each subdomain Ωi
196
is generated using a suitable meshing algorithm and the most appropriate element type. The resulting meshes should
197
not share nodes, edges, faces, or cells. Finally, we combine these submeshes to obtain the non-conformal hybrid
198
mesh of the whole domain Ω. This strategy of non-conforming hybrid mesh generation has several advantages over
199
other approaches based on hybrid conforming mesh. For instance, it allows one (i) to avoid the intermediate layer
200
of pyramids as transition elements between triangular and quadrilateral faces; (ii) to relax the node-to-node matching
201
requirement at the interface between the meshes; (iii) to change or replace certain portions of the domain without
202
changing the entire mesh. These features are particularly useful for parametric modelling of complex geometries
203
where the object is built from several parts designed separately. Figure 2(a)-(b) illustrates examples of hybrid meshes
204
generated by this strategy.
205
After combining the meshes, we need to find the hybrid non-conforming interfaces between triangular and quadri-
206
lateral faces in order to compute the flux matrices, see Figure 2(c). First, we locate the boundary faces of each
207
subdomain mesh and store them in a list including their element and face location. Next, we use the fact that only 10
Page 10 of 28
the boundary faces of two subdomains that shared boundaries may participate in a partial face-face connection. Let
209
us assume that two subdomains Ωi and Ω j share a boundary. For each of the boundary faces of the mesh of Ωi , we
210
perform two tests: the first is an intersection test with every member of the boundary face list of the mesh of Ω j . This
211
just indicates whether or not an intersection exists. The second is a find-intersection test that involves computing the
212
set of intersections between triangular and quadrilateral faces. To perform these intersection tests, we use the class
213
Intersection provided by the Eberly’s Wild Magic library [22]. For the reader’s convenience, we briefly describe
214
the general idea of this class. The test-intersection is based on the method of separating axes [34], a method for deter-
215
mining whether or not two convex objects are intersecting. A test for nonintersection of two convex objects is simply
216
stated: if there exists a line for which the intervals of projection of the two objects onto that line do not intersect, then
217
the objects do not intersect. Such a line is called a separating axis. For a pair of convex polygons in 3D, only a finite
218
set of direction vectors needs to be considered for separation tests. The polygons normals are potential separating
219
directions. If either normal direction separates the polygons, then no intersection occurs. However, if neither normal
220
direction separates the polygons, then two possibilities exist for the remaining potential separating directions. The
221
first case is that the polygons are coplanar and the remaining potential separating directions are those vectors in the
222
common plane and perpendicular to the polygon edges. The second case is that the polygon planes are not parallel
223
and that each polygon is split by the plane of the other polygon. The remaining potential separating directions are
224
the cross products of two edges, one from each polygon. In summary, for the intersection between a triangular and
225
a quadrilateral face in 3D, a maximum of 19 separating directions need to be tested in the worst case. Given a line,
226
it is straightforward to project the vertices of a polygon onto the line and compute a bounding interval for those pro-
227
jections. The two intervals obtained from the two polygons are easily compared for overlap. If an intersection exist,
228
we perform a find-intersection query using Boolean operations on polygons [47, 39]. Here the intersection between
229
a triangular and a quadrilateral face (both lies in the same plane) is a convex polygon with at most seven vertices.
230
Let Q be a quadrilateral and T be a triangle with vertices (v1 , v2 , v3 ) that are counterclockwise ordered. We define
231
half-planes, πi j , 1 ≤ i < j ≤ 3, bounded by the lines (vi v j ) and containing the third vertex of the triangle T . The
232
intersection between T and Q is defined as T ∩ Q = ((Q ∩ π12 ) ∩ π23 ) ∩ π13 . The problem is reduced to find the
233
intersection between a polygon and a half-plane which is relatively straightforward. For a more complete exposition
234
of the theory of intersection tests, we refer the reader to the books [49, Chap. 7 & 11]-[21, Chap. 15].
235
4.2. Non-conforming hybrid interface matrix calculation
Ac ce pt e
d
M
an
us
cr
ip t
208
236
A key ingredient for the implementation of the DG method on hybrid meshes is the computation of the matrix
237
Aτh involving integrals over a hybrid interface. We concentrate here on the evaluation of such integrals. In fact, the
238
other matrices in Eq. (9) which involve integrals over an element or an interface between two elements of the same
239
type are evaluated as specified in Eq. (5). We just have to utilize the correct set of basis functions corresponding to
240
the particular element type.
241
Let aik = τi ∩ hk be a hybrid interface between a tetrahedron τi and a hexahedron hk . Here aik is a n-gon with 11
Page 11 of 28
1111111 0000000 ■ 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111
.
Hy
.
fo on −c on es d n ac bri interf
■
ip t
ng
(c)
cr
(b)
..
.
i rm
(a)
111111 000000 000000 111111 000000 111111 000000 111111 000000 111111 00000000000000 11111111111111 000000 111111 00000000000000 11111111111111 000000 111111 00000000000000 11111111111111 000000 111111 00000000000000 11111111111111 000000 111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 ■11111111111111 00000000000000 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 ■
Figure 2: Examples of non-conformal hybrid meshes generated using the combining meshes strategy. In (a) we show a triangular-quadrilateral
us
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
an
node (filled circle).
242
n = 3, 4, 5, 6, 7, that is, aik is a convex polygon with n edges (see Figure 3). Let also R be a regular polygon in R2 and
243
Υ : R ∋ ξ → x ∈ aik ⊂ R3 be an explicit parametrization of aik . The matrix Aτh is block diagonal, whose non-zero blocks are equal to the interface matrix Sik when aik is a hybrid interface. The matrix Sik is defined as Z Z ) −1 k i) ξ ) L(q ξ ) k∇Υ(ξξ )kdξξ L(p ◦ Ψ−1 θi j (x)ϑkl (x)dx = (Sik ) jl = hk ◦ Υ(ξ j ◦ Ψτi ◦ Υ(ξ l R
aik
M
244
(13)
where θi j ∈ P pi (τi ) and ϑkl ∈ Qqk (hk ) are the basis functions of degrees pi and qk in the elements τi and hk , respec-
246
tively. The integrals in Eq. (13) are calculated numerically and stored for each hybrid interface aik . To this end, two
247
approaches are possible
Ac ce pt e
d
245
248
1. Partitioning the interface aik (a n-gon) into n triangles and then performing cubature rules over each triangle
249
[55, 17]. For n-gons with n ≥ 4, we use the centroid of aik to partition it into triangles. The matrix Sik is
250
computed with the following formula
(Sik ) jl =
251
252
253
254
θi j (x)ϑkl (x)dx =
aik
s=n Z X s=1
θi j (x)ϑkl (x)dx,
(14)
∆s
where ∆ s is a triangle formed by two successive vertices of aik and its centroid as illustrated in Figure 4(a). The P parametrization of ∆ s can be defined as x = Υ(ξξ ) = 3i=1 φi (ξξ )vi , where ξ = (ξ1 , ξ2 ) ∈ T 2 , φi are the barycentric coordinates in T 2 and vi denotes the three vertices of ∆ s . To compute the last integral in Eq. (14), we lay down Nc cubature points λn ∈ T 2 and weights ωn . This is achieved by Z
∆s
255
Z
θi j (x)ϑkl (x)dx =
Nc X
(p )
(q )
L j i (λτn )Ll k (λhn )ωn k∇Υ(λn )k,
(15)
n=1
3 h −1 3 where λτn = Ψ−1 τi ◦ Υ(λn ) ∈ T and λn = Ψhk ◦ Υ(λn ) ∈ D .
12
Page 12 of 28
2. An alternative is to use cubature rules over polygons. However, the design and development of such rules has
257
not yet reached a mature stage. To our knowledge this approach has not been used previously in the context
258
of DG methods. Recently, Mousavi et al. [43] presented an algorithm to build symmetric cubature rules for
259
arbitrary regular polygons. This algorithm combines elements of group theory and numerical optimization and
260
was originally proposed by Wandzura and Xiao [55] for the construction of cubature rules over triangles. Let Pn
261
be a reference regular n-gon with vertices pm = (cos 2πm/n, sin 2πm/n) for m = 1, 2, . . . , n in the ξ -coordinate
264
point p ∈ Pn defined as [53]
αi (ξξ ) , φi (ξξ ) = Pn ξ) j=1 α j (ξ
α j (ξξ ) =
(1 − |ξξ |2 ) sin3 (π/n) cos(π/n) , Ai−1 Ai
cr
263
system, see Figure 4(b). The parametrization Υ from Pn to aik is obtained via the isoparametric mapping P x = Υ(ξξ ) = ni=1 φi (ξξ )vi , where vi are the vertices of the n-gon aik and φi are the Laplace shape functions at ξ = (ξ1 , ξ2 ) ∈ Pn ,
us
262
ip t
256
where αi (ξξ ) is the Laplace weight function, Ai is the area of the triangle [p, pi , pi+1 ] and Ai−1 is the area of the
266
triangle [p, pi−1 , pi ]. The Laplace shape functions are non-negative, interpolate nodal data and form a partition
267
of unity. They are also linear on the polygon boundaries which ensures that the parametrization Υ(ξξ ) is an affine
268
mapping. Once Υ(ξξ ) is given explicitly, the integral in Eq. (13) is evaluated by laying down Nc cubature points
269
on Pn and proceed as in Eq. (15).
M
an
265
We compare the two approaches by evaluating the integral of x6 y6 z6 over the regular pentagon, hexagon and heptagon.
271
The relative errors between exact and approximate values are shown in Figure 5 as a function of the number of cubature
272
points. We observe that for a relative error of 10−8 the polygonal cubatures require 3 to 5 times fewer points than the
273
partitioning approach. Hence, for the numerical results given in Section 5, we adopt the second approach since it is
274
more practical and requires the least number of cubature points for a given degree of accuracy.
Ac ce pt e
d
270
●
●
●
aik
aik
●
●
●
●
aik
●
●
●
● ● ●
● ●
●
●
aik
● ●
●
●
aik
● ●
● ●
Figure 3: Five possible configurations for the hybrid interface aik between a quadrilateral face and a triangular face in 3D.
275
5. Numerical results
276
5.1. Two-dimensional examples
277
We consider here the 2D TMz -polarized Maxwell equations, that is, we solve Eq. (1) for the fields H = (H x , Hy , 0)
278
and E = (0, 0, Ez). Two examples with analytical solutions are studied in detail to illustrate the accuracy and per-
279
formance of the DG method on hybrid triangular-quadrilateral meshes. The exact time-domain solutions of these
280
problems can be found in [24]. 13
Page 13 of 28
. . ξ2
p2
.
T2
ξ1
.
v3
x1
.
3
v2
∆s
.
■
.
x3
.
p
.
x2
1
.
p6
ξ1
.
. . . . P6
p4
.
v1
p
p
v2
v3
5
x2
x1
(a)
v1
aik
x3
v4
.
ip t
ξ2
.
.
v6
cr
. .
v5
(b)
us
Figure 4: Two numerical integration schemes used for the evaluation of integrals over a polygon aik . In (a) we show a scheme based on the partition of aik into triangles ∆ s . The scheme illustrated in (b) consists of using cubature rules over regular polygons. 100 Hexagonal cubature Partitioning approach
10-4 135 points
198 points
10-8 47 points 0
50
100 150 200 Number of cubature points
47 points
46 points 250 0
50
100 150 200 Number of cubature points
d
10-10
231 points
M
Relative error
10-2
10-6
Heptagonal cubature Partitioning approach
an
Pentagonal cubature Partitioning approach
250 0
50
100 150 200 Number of cubature points
250
Ac ce pt e
Figure 5: Convergence curves of two numerical integration approaches used to evaluate the integral of x6 y6 z6 over the regular pentagon (left), hexagon (center) and heptagon (right).
281
5.1.1. Eigenmode in a PEC cavity
282
In the first example, we propose to validate the implementation of the proposed DGTD-P p Qq method on non-
283
conforming hybrid meshes by performing a numerical convergence test. To this end, we consider the propagation of
284
an eigenmode which is a standing wave of frequency f = 212 MHz and wavelength of 1.4 m in a unitary vacuum-filled
285
square cavity with PEC walls. The computations are performed on three successively refined hybrid meshes, named
286
MSi for i = 1, 2, 3. Each mesh MSi consists of 5 × 22i−1 triangles and 3 × 22i−1 quadrilaterals, as shown in Figure 6.
287
We present results with interpolation orders, p and q, up to degree 8 such that p = q. For each interpolation order,
288
we use the maximum time step, ∆t, allowed by the stability condition of the hybrid method. In Figure 7 (left) we
289
illustrate the h-convergence graphs as a function of the total number of degrees of freedom (#DOF). As expected, for
290
smooth solutions the convergence rate behaves as O(∆t4 + h p ) even for non-conforming hybrid meshes. Also shown in
291
Figure 7 (right) is the p-convergence rate in linear-log scale, observing an exponential decay of the error as exp(−s p p),
292
with decay constant s p = log[E(p − 1, q − 1)/E(p, q)] ≃ 2, where E(p, q) is the L2 error when the interpolation orders
293
p and q are used. By inspecting the results in Figure 7, we observe that the gain in accuracy is notable when refining 14
Page 14 of 28
the mesh or increasing the interpolation orders. Furthermore, to achieve a certain accuracy level, the DG methods
295
require less #DOF for higher interpolation orders, as typically found with other DG methods based on fully simplicial
296
meshes [25, 30]. 1
0.75
0.75
0.75
0.5
0.5
0.5
0.25
0.25
0.25
0
0 0
0.25
0.5
0.75
cr
1
us
1
ip t
294
0
1
0
0.25
0.5
0.75
1
0
0.25
0.5
0.75
1
an
Figure 6: Sequences of three successively refined hybrid meshes used for assessing convergence of the eigenmode problem.
100
10-2
10-2
p=q=2, rate=2.08 p=q=3, rate=2.92
10-6
L2 error
10-4
10-4
p=q=4, rate=3.87
d
L2 error
M
p=q=1, rate=1.25
p=q=5, rate=4.10
10-8
10-6
Ac ce pt e
p=q=6, rate=4.24
p=q=7, rate=4.08
p=q=8, rate=4.13
10-10
10
10-8 2
100
(DOF)1/2
3
4
5 6 Interpolation order
7
8
Figure 7: On the left is shown the L2 -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(∆t4 + h p ) for all orders. On the right, we plot the slope of the error for the mesh MS1 , confirming the spectral convergence of the hybrid DGTD-P p Qq method.
297
5.1.2. Plane-wave scattering by a dielectric circular cylinder
298
As an illustration of the capability of the hybrid DG method to handle materials, let us consider plane wave
299
scattering by a penetrable circular cylinder with a radius of 0.6 m consisting of a dielectric nonmagnetic material
300
with a relative permittivity εr = 12.0. The frequency of the incident wave is f = 300 MHz and the wavelength is
301
λ = 0.288 m. The surrounding medium is assumed to be vacuum. The computational domain Ω is bounded by a
302
square of side length 3.2 m centered at (0, 0). A Silver-Müller absorbing condition is applied on the boundary of
303
the square. We compare the performance of the DGTD-P p Qq method using a non-conforming hybrid mesh with the
304
DGTD-P p method on a fully triangular mesh, see Figure 8. The hybrid mesh is generated by combining submeshes 15
Page 15 of 28
of three subdomains, Ω1 = [−0.35, 0.35]2 m2 (inside the cylinder), Ω2 = [−0.7, 0.7]2 m2 \Ω1 (close to cylinder and
306
completely enclosing it) and Ω3 = [−1.6, 1.6]2 m2 \(Ω1 ∪ Ω2 ) such that Ω = Ω1 ∪ Ω2 ∪ Ω3 . The regions Ω1 and
307
Ω2 are filled with regular quadrilaterals while Ω3 is filled with unstructured triangles. The global hybrid mesh has
308
1888 triangles, 1204 quadrilaterals, 5044 conforming interfaces, 320 non-conforming hybrid interfaces and 120 non-
309
conforming interfaces between quadrilaterals. The average edge length of this mesh is of 0.16λ. The triangular mesh is
310
obtained by replacing the quadrilateral elements of the hybrid mesh by unstructured triangles. The resulting triangular
311
mesh is conformal and consists of 4376 triangles and 6626 interfaces with an average edge length of 0.18λ. For the
312
hybrid DGTD-P p Qq method, the interpolation order, q, is a pair q = (q1 , q2 ), where q1 and q2 are the interpolation
313
orders used in the smallest and biggest quadrilaterals. The final simulation time has been set to T = 33.33 ns which
314
corresponds to 10 periods of the incident wave. The results are listed in Table 1 in terms of accuracy and computational
315
cost to reach the final simulation time for different p and q. Inspecting these results, we observe that for a fixed p,
316
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
317
the triangular method. On Figure 9 we compare the x-wise 1D distributions of the real part of the discrete Fourier
318
transform of Ez and Hy components. We observe from this figure that, the main error of DG solutions arises at the
319
material interfaces where the solution loses smoothness and the results confirm the ability of the hybrid DG method to
320
accurately model problems with discontinuous coefficients and solutions. Finally, contour plots of the Ez component
321
at time t = T for numerical simulations performed with the DG method on the hybrid and the triangular meshes, are
322
illustrated in Figure 10 for p = 2 and q = (2, 2). We observe very good agreement between the solutions obtained
323
by both DG methods and the exact solution even with relatively low grid resolutions of six points per λ. The hybrid
324
method shows slightly better results than the triangular method.
Ac ce pt e
d
M
an
us
cr
ip t
305
1.6
1.6
0.8
0.8
0
0
-0.8
-0.8
-1.6
-1.6 -1.6
-0.8
0
0.8
1.6
-1.6
-0.8
0
0.8
1.6
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.
16
Page 16 of 28
Table 1: Comparison between the DGTD-P p Qq 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 (q1 , q2 ) where q1 and q2 are the interpolation orders used in the smallest and biggest quadrilaterals, respectively. We give the maximum ∆t allowed by the stability condition, the L2 -errors at time t = 33.33 ns, the #DOF and the
L2 -error
#DOF
CPU (s)
DGTD-P2
3.648E-02
2.394E-01
26256
201
DGTD-P3
2.132E-02
3.931E-02
43760
581
DGTD-P4
1.684E-02
6.532E-03
65640
1108
DGTD-P5
1.403E-02
1.083E-03
91896
1886
DGTD-P6
1.038E-02
1.785E-04
122528
3488
DGTD-P2 Q(2,2)
3.648E-02
1.692E-01
22164
165
DGTD-P3 Q(3,2)
2.346E-02
2.398E-02
32768
379
DGTD-P4 Q(3,2)
1.841E-02
5.105E-03
42208
592
DGTD-P5 Q(4,3)
1.523E-02
5.784E-04
62836
1162
DGTD-P6 Q(4,3)
1.216E-02
1.484E-04
76052
1671
cr
an
Method
ip t
∆t (ns)
us
computing time (CPU) in seconds to achieve the final time period.
Exact Ez DG-P3 DG-P3Q(3,2)
1
M
8
1.5
Exact Hy DG-P3 DG-P3Q(3,2)
6 4
-0.5
-1
-1.5 -1.8
d -0.6
x
0.6
field Hy
0
Ac ce pt e
field Ez
0.5
2 0 -2 0.7
zoom-in
1.2
zoom-in
-4 -6 -0.1 -0.7
-8 -1.8
1.8
0.9 0.6
-0.6
-0.6
x
0.6
0.7
1.8
Figure 9: Plots of the Ez (left) and Hy (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-P3 Q(3,2) (square) and triangular DGTD-P3 (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.
325
5.1.3. Standing wave in a wedge-shaped PEC resonator
326
As an example with geometric singularities, we consider the propagation of a standing wave of frequency f =
327
1293 MHz and wavelength λ = 0.232 m inside a domain bounded by the curves y = tan(7π/4)x, x2 + y2 = 1/4 and the
328
x-axis. The non-convex nature of the outer boundary renders some components of the solution highly singular, making
329
it a challenging problem. We compare approximate solutions obtained using the DGTD-P p method formulated on 17
Page 17 of 28
1
0.5
0.5
0.5
0
0
0
-0.5
-0.5
-0.5
-1 Y -1 Z X
-0.5 0 0.5 Ez component - exact solution
1
-1 Y -1 Z X
-0.5 0 0.5 Ez component on hybrid mesh for p=2, q=(2,2)
1
-1 Y -1 Z X
ip t
1
-0.5 0 0.5 Ez component on triangular mesh for p=2
cr
1
1
Figure 10: Contour lines of the Ez component at time t = 33.33 ns for the exact solution (left) and for solutions resulting from the hybrid DGTD-
us
P2 Q(2,2) (center) and triangular DGTD-P2 (right) methods. Zoomed-in region to show more details. In all the plots, we use 13 equally spaced contours from Ez = −1.5166 to Ez = 1.5299.
triangular meshes with those resulting from the DGTD-P p Qq method formulated on hybrid non-conforming meshes.
331
The triangular mesh consists of 6280 elements and 9531 interfaces with an average edge length of 0.057λ. The hybrid
332
mesh contains 1745 triangles, 1820 quadrilaterals, 6228 conforming interfaces and 280 non-conforming interfaces.
333
The global hybrid mesh has an average edge length of 0.082λ. Moreover, we can see that the hybrid mesh is locally
334
refined in the vicinity of the re-entrant corner located at the origin. The resulting meshes are depicted in Figure 11.
335
The simulation time has been set to T = 38.5 ns which corresponds to 50 periods. Results are summarized in Table 2
336
in terms of accuracy and computational costs. Here, the hybrid DG method produces numerical solutions that are one
337
order of magnitude more accurate than those delivered by the triangular DG method. Figure 12 shows the contour lines
338
of the H x component at t = T for numerical solutions obtained using the DGTD-P4 and DGTD-P2 Q4 method. One
339
can clearly observe a good qualitative agreement between the exact solution and the hybrid DGTD-P2 Q4 solutions.
340
The triangular DGTD-P4 solution produces some oscillations near the re-entrant corner due to bad resolution. This
341
is confirmed in Figure 13 where we plot the x-wise 1D distributions of the H x and Hy components along the line
342
y = 5x + 0.025 passing nearby the corner singularity. Inspecting Table 2 and Figure 13, we observe that, for the
343
DGTD-P p method, the error decreases slowly when increasing the order from p = 2 to p = 4. In contrast, the hybrid
344
method offers faster convergence when increasing the order in the quadrilateral elements from q = 2 to q = 4. More
345
specifically, the global error depends on the global regularity of the exact solution which motivates the use of lower
346
polynomial order together with mesh refinement near singular points in order to recover standard convergence rates.
347
5.2. Three-dimensional examples
348
5.2.1. Eigenmode in a spherical cavity
Ac ce pt e
d
M
an
330
349
As a first verification of the general 3D framework, let us consider the propagation of the (0, 1, 1)-mode inside
350
a PEC spherical cavity of unit radius. The resonant frequency is 131 MHz and the wavelength is λ = 2.29 m. The
351
exact time-domain solution of this problem can be found in [3]. We compare the accuracy and efficiency of the DG 18
Page 18 of 28
0.5
0
0
-0.5 -0.5
-0.5 -0.5
0.5
0
0.5
us
0
cr
ip t
0.5
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.
an
Table 2: Comparison between the hybrid DGTD-P p Qq method and the triangular DGTD-P p method. We give the maximum ∆t allowed by the
L2 -error
#DOF
CPU (s)
DGTD-P2
1.570E-02
5.482E-02
37680
1604
DGTD-P4
7.246E-03
2.381E-02
94200
9832
DGTD-P2 Q2
1.724E-02
9.875E-03
26850
1211
DGTD-P2 Q4
8.701E-03
1.043E-03
55970
6539
Method
-0.18 -0.18
0.18
Ac ce pt e
0.18
d
∆t (ns)
M
stability condition, the L2 -errors at time t = 38.5 ns, the #DOF and the computing time (CPU) in seconds to achieve the final time period.
0.18
Hx component - exact solution
-0.18 -0.18
0.18
0.18
Hx component on hybrid mesh for p=2, q=4
-0.18 -0.18
0.18
Hx component on triangular mesh for p=4
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 DGTDP2 Q4 (center) and triangular DGTD-P4 (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.
352
method using hexahedral, tetrahedral and hybrid meshes. First, a fully hexahedral mesh consisting of 6413 hexahedra
353
and 19602 quadrilateral faces is generated. Then, the fully tetrahedral mesh is obtained by considering the hexahedral
354
mesh and splitting each hexahedron to six tetrahedra. The resulting tetrahedral mesh has 38478 tetrahedra and 77682 19
Page 19 of 28
0.25
0.4
Exact Hx DGTD-P2Q4 DGTD-P4
0.3
Exact Hy DGTD-P2Q4 DGTD-P4
0.2
0.2
zoom-in 0.11
0.15 0.1 -0.004
field Hy
field Hx
0.1 0
0
0.1 0.05
ip t
-0.1 0
-0.2
-0.05
-0.3 -0.4
-0.1 -0.05
0
0.05
0.1
-0.1
-0.05
0
0.05
0.1
cr
-0.1
x
x
Figure 13: Plots of the H x (left) and Hy (right) components along the line y = 5x + 0.025 at time t = 38.5 ns for the wedge-shaped cavity. We
us
compare the exact solution (solid line) with solutions obtained using the DGTD-P2 Q4 (dashed line) and DGTD-P4 (dotted line) methods.
triangular faces. Finally, the hybrid mesh contains 13068 tetrahedra, 216 regular hexahedra, 26136 triangular faces,
356
540 quadrilateral faces and 2916 hybrid non-conforming interfaces. These meshes are illustrated in Figure 14. The
357
average edge lengths of the hexahedral, tetrahedral and hybrid meshes are respectively 0.0387λ, 0.0463λ and 0.0676λ
358
which correspond to 26, 22 and 15 points per wavelength. We report on results obtained by the hexahedral DGTD-Qq ,
359
tetrahedral DGTD-P p and hybrid DGTD-P p Qq methods with interpolation orders p and q ranging from p = 1 up to
360
p = 5 with p = q. The simulation time is fixed to T = 76 ns which corresponds to 10 periods. Table 3 summarizes the
361
values of time steps used in the simulation. In Figure 15 (left) we show the L2 error as a function of the interpolation
362
order. Inspecting these results we observe several things. First, the gain in accuracy is notable when increasing the
363
interpolation order for all methods. Here, the exponential decay of the error as exp(−2p) is observed for all cases.
364
Second, the tetrahedral method is 3.9 to 7.9 times more accurate than the hexahedral method when increasing the
365
interpolation order from p = 1 to p = 5. The results reveal that the accuracy of the hybrid method is almost identical
366
to that of the tetrahedral method. Figure 15 (center) shows the ratio between the #DOF required by the hexahedral,
367
tetrahedral and hybrid methods to achieve the final simulation time. The results indicate that the hybrid method is able
368
to reduce the #DOF for given p and q. For instance, with p = q = 5, the L2 error is 7.1 × 10−5 for the hybrid mesh,
369
4.7 × 10−4 for the hexahedral method and 6.0 × 10−5 for the tetrahedral mesh. With the hybrid mesh, the #DOF is
370
reduced by a factor of 1.8 compared to hexahedral mesh and by factor of 2.8 compared to tetrahedral mesh. Note that,
371
the higher #DOF required by the tetrahedral mesh stems from the fact that the number of elements in the tetrahedral
372
mesh is 6 times higher than those in the hexahedral mesh. Finally, the comparison in terms of computing times is
373
displayed in Figure 15 (right). It can be observed that the reduction in the #DOF is translated into a corresponding
374
reduction in computational cost. Here, the hybrid method allows for a reduction in computing time by a factor ranging
375
from 2.4 to 8.0 compared to hexahedral method and by a factor of 3.5 compared to tetrahedral method. The lower
376
computing times of the hybrid method partially stems from its larger time steps employed in hybrid mesh compared
377
to those used in tetrahedral and hexahedral meshes.
Ac ce pt e
d
M
an
355
20
Page 20 of 28
ip t cr
Y X
Y Z
X
Y
Z
Z
X
us
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.
an
Table 3: The time steps allowed by the stability conditions of the DG method using hexahedral, tetrahedral and hybrid meshes for different
p=q=2
Hexahedra
1.970E-02
9.850E-03
Tetrahedra
2.357E-02
1.178E-02
Hybrid
3.178E-02
1.589E-02 3
HEX TET HYB
Ratio of the DOF
10-2 10-3
10
-4
10
-5
p=q=5
5.910E-03
4.728E-03
3.546E-03
7.072E-03
5.658E-03
4.243E-03
9.535E-03
7.628E-03
5.721E-03
8
HEX/TET HEX/HYB
Ac ce pt e
L2 error
10-1
p=q=4
2
Ratio of computing times
100
p=q=3
d
p=q=1
M
interpolation orders.
TET/HYB
1
2
3
Interpolation order
4
5
1
HEX/TET HEX/HYB
6
TET/HYB
5 4 3 2 1
0
1
7
0
2
3
4
Interpolation order
5
1
2
3
4
5
Interpolation order
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.
378
5.2.2. Exposure of human head tissues to a localized source radiation
379
As a final example, we consider a more challenging problem which is concerned with the propagation of an EM
380
wave emitted by a localized source in a realistic geometrical model of head tissues. Starting from magnetic resonance
381
images of the Visible Human 2.0 project [46], head tissues are segmented and the interfaces of a selected number of
382
tissues are triangulated. Different strategies can be used in order to obtain a smooth and accurate segmentation of 21
Page 21 of 28
head tissues and interface triangulations as well. The strategy adopted in this example consists in using a variant of
384
Chew’s algorithm [9], based on Delaunay triangulation restricted to the interface, which allows to control the size and
385
aspect ratio of interface triangles. Example of triangulations of the skin, skull and brain are shown in Figure 16. Then,
386
these triangulated surfaces are used as inputs for the generation of volume meshes of the tissues. Finally, the GHS3D
387
tetrahedral mesh generator [33] is used to mesh the volume domains between the various interfaces. The exterior of
388
the head must also be meshed, up to a certain distance from the skin. The computational domain is here artificially
389
bounded by a surface, denoted by Γa , on which the absorbing boundary condition is imposed. Overall, the constructed
390
geometrical model considered here consists of four tissues (skin, skull, CSF - Cerebro Spinal Fluid and brain). The
391
characteristics of the tissues are summarized in Table 4 where the values of the relative permittivity εr , and the
392
conductivity σ, correspond to a frequency F = 1800 MHz and have been obtained from a special purpose online data
393
base. For all tissues, the relative permeability, µr , is set to 1. Given this tetrahedral mesh of the head tissues, the hybrid
394
tetrahedral-hexahedral mesh is built in the following way (see the left panel of Figure 17): (a) the artificial boundary
395
Γa is defined as a plane-parallel surface; (b) a second plane-parallel surface, denoted by Γ0 , lying close to the head
396
and completely enclosing it, is introduced; (c) the volume space between Γa and Γ0 (vacuum space) is discretized
397
using a regular hexahedral mesh; (d) the volume space between Γ0 and the skin is discretized with unstructured
398
tetrahedral elements. The global hybrid mesh consists of 197691 tetrahedra, 7000 hexahedra, 394397 triangular
399
interfaces, 21900 quadrilateral interfaces and 5124 hybrid non-conforming interfaces. The minimum, maximum and
400
average edge lengths of the hybrid mesh are equal to 1.12 mm, 43.14 mm and 17.26 mm, respectively. To compare the
401
performance of the DG methods on hybrid and tetrahedral meshes, a fully tetrahedral mesh is generated by considering
402
the hybrid mesh and replacing the hexahedral part by unstructured tetrahedral elements. We stress that the resulting
403
tetrahedral mesh is conforming i.e., the same triangulated surface Γ0 is used for the generation of tetrahedral elements
404
inside and outside Γ0 , as shown in Figure 17 (right). The global tetrahedral mesh consists of 350581 tetrahedra and
405
709395 interfaces. The minimum, maximum and average edge lengths of the tetrahedral mesh are equal to 0.49 mm,
406
44.32 mm and 13.98 mm, respectively. Finally, a dipolar type source is localized near the right ear of the head
407
yielding a current of the form Jz (x, t) = δ(x − xd ) f (t), where f (t) is a sinusoidally varying temporal signal and xd is
408
the localization point of the source. The physical simulation time has been fixed to T = 2.78 ns which corresponds to
409
five periods of the temporal signal. A discrete Fourier transform of the components of the electric field is computed
410
during the last period of the simulation. For this problem, we consider using the DGTD-P2 method on the tetrahedral
411
mesh and the DGTD-P2 Q2 and DGTD-P2 Q1 methods on the hybrid mesh. The time steps used in the simulations
412
are ∆t = 2.728E-04 ns for the DGTD-P2 method, ∆t = 6.655E-04 ns for the DGTD-P2 Q2 method and ∆t = 1.012E-
413
03 ns for the DGTD-P2 Q1 method. These values correspond to the maximum ∆t allowed by the stability conditions
414
of the DG methods. Contour lines of the modulus of the electric field, |E|, on selected planes of the skin for the
415
approximate solutions resulting from the hybrid DGTD-P2 Q2 and tetrahedral DGTD-P2 methods are visualized in
416
Figure 18. Visually both methods provide almost the same solutions except slight amplitude differences around the
417
right ear of the head. For the sake of completeness, we compare in Figure 19 the time evolution of the Ez and Hy
Ac ce pt e
d
M
an
us
cr
ip t
383
22
Page 22 of 28
components at a selected point in the center of the head. The comparisons show a visually perfect match between all
419
solutions and no numerical artifacts due to the non-conformity of the mesh is observed. In Figure 19 we also plot
420
the differences between the DGTD-P2 solution and the hybrid solutions, showing that the errors are relatively very
421
small (less than 0.1% for DGTD-P2 Q1 and less than 0.01% for DGTD-P2 Q2 ). Finally, we give some informations on
422
the simulation times required by both schemes. On a workstation equipped with an Intel Xeon CPU E5-2650 @2.00
423
GHz and 64 GB of RAM, the DGTD-P2 method requires 27 h 47 min for a total of 10183 time steps, the DGTD-P2 Q2
424
method requires 7 h 35 min for a total of 4177 time steps, and the DGTD-P2 Q1 method requires 4 h 17 min for a total
425
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
426
DG method.
M
an
us
cr
ip t
418
Z
Z
YX
Z Y
X
d
YX
Ac ce pt e
Figure 16: Surface mesh of the skin (left), the skull (center) and the brain (right).
Table 4: Electromagnetic characteristics of the selected head tissues at frequency f = 1800 MHz.
Tissue Skin Skull CSF Brain
427
εr
σ (S/m)
wavelength (mm)
38.87
1.18
26.73
15.56
0.43
42.25
67.20
2.92
20.33
43.55
1.15
25.26
6. Concluding remarks
428
We have reported on our efforts to design a high-order accurate DG method for solving the time-domain Maxwell
429
equations on non-conforming hybrid meshes. The proposed method combines a central numerical flux with a fourth-
430
order leap-frog time integration. Within each mesh element, the EM field components are approximated using ar-
431
bitrary high-order nodal polynomials. The geometrical flexibility is achieved by using unstructured triangular or 23
Page 23 of 28
ip t cr us
an
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.
tetrahedral meshes around complex geometric details and computational efficiency is then improved by using quadri-
433
laterals or hexahedra to fill the remainder of the domain. The resulting DGTD method conserves a discrete form of
434
the EM energy, and is stable under some CFL-like conditions. The key issues for the practical implementation of
435
the proposed method are the generation of the hybrid meshes and the computation of flux matrices at the interfaces
436
between hexahedral and tetrahedral faces. First, the mesh strategy adopted in this paper relies on the decomposition
437
of the computational domain into non-overlapping subdomains, each one is meshed independently using the appro-
438
priate element type. Then, the hybrid interfaces are found using the method of separating axis and Boolean operations
439
over polygons. This strategy demonstrates the benefits of relaxing the conformity requirement in the discretization
440
process in view of facilitating the construction of meshes for complex propagation configurations, and improving the
441
overall efficiency of the simulations. It is noteworthy that, most, if not all, of commercial or non-commercial mesh-
442
generating softwares currently produce only conforming hybrid meshes by adding transition layers of pyramids or
443
prisms. Second, the flux matrices are calculated using cubature rules over polygons. For a given accuracy, these rules
444
are able to save more than 65% of the cubature points compared to other approaches. The proposed DGTD method is
445
implemented in two and three space dimensions and tests on different wave propagation problems involving complex
446
realistic geometries and heterogeneous media are presented. The spectral convergence for p-refinement method and
447
the optimal convergence rate for h-refinement method have been verified. Results also show excellent agreement with
448
exact solutions and no numerical artifacts caused by the non-conformity of the mesh are observed. Compared to DG
449
methods on single element-type meshes, the DG method on hybrid meshes leads to clear reduction in computing time
450
which is highly dependent on the interpolation order. In 2D, it is found that the hybrid mesh is able to produce at least
451
a similar accuracy and allows the DG method to save up to 55% of the computing time and the number of degrees of
Ac ce pt e
d
M
432
24
Page 24 of 28
0.514
0.481
0.481
0.447
0.447
0.413
0.413
0.38
0.38
0.346
0.346
0.313
0.313
0.279
0.279
0.245
0.245
0.212
0.212
ip t
0.514
0.178
0.178
0.144
0.144
0.111
0.111
0.0772
0.0772
0.0436
0.504 0.471 0.438 0.406 0.373 0.307 0.274 0.241 0.209
0.0772 0.0443 0.0115
cr
Y Z X
|E| - DGTD-P2Q2
0.504 0.471 0.438 0.406 0.373 0.34 0.307 0.274 0.241 0.209 0.176 0.143 0.11 0.0772 0.0443 0.0115 |E| - DGTD-P2
d
Y Z X
M
0.176 0.11
|E| - DGTD-P2
an
0.34
0.143
0.01
Z Y X
|E| - DGTD-P2Q2
us
Z Y X
0.0436
0.01
Ac ce pt e
Figure 18: Contour lines of the modulus of the electric field, |E|, in selected cut planes of the head, using the hybrid DGTD-P2 Q2 method (left) and the tetrahedral DGTD-P2 method (right).
452
freedom. In 3D, the gain becomes more significant, as the computational cost is reduced by up to 85%. A comparison
453
between a fully hexahedral mesh and a fully tetrahedral mesh (obtained by converting each hexahedron to six tetrahe-
454
dra) reveals that the tetrahedral method is always more accurate and more efficient than hexahedral method, excepting
455
for linear interpolation orders. The parallelization of the proposed hybrid DG method and the comparison with DG
456
methods on conforming hybrid meshes are left for future study.
457
Acknowledgements
458
The author would like to thank Dr. Seyed Mousavi of the University of California at Davis, USA, for providing
459
the quadrature points over polygons to calculate flux matrices at hybrid interfaces. We also thank Dr. Nicolas Roquet
460
for several discussions concerning the intersection tests for hybrid mesh generation. Finally, the author acknowledge
461
partial support from the French National Agency for Research (ANR-10-JCJC-0905).
25
Page 25 of 28
0.05
0.1 zoom-in
0.043
0.04
0.08
0.03
0.06 0.04
0.01 2.57
2.65
0 -0.01
0.02 0 -0.02
-0.02
DGTD-P2
-0.03
DGTD-P2Q2
-0.04
DGTD-P2
-0.06
DGTD-P2Q2
DGTD-P2Q1
DGTD-P2Q1 -0.08
cr
-0.04
-4
-3
10
10
(P2 - P2Q1)
0
0
(P2 - P2Q2) 0
0.5
(P2 - P2Q2)
1 1.5 2 time (nanoseconds)
-10
2.5
-4
0
0.5
an
-3
us
error
error
(P2 - P2Q1)
-10
ip t
0.037
field Ez
field Hy
0.02
1 1.5 2 time (nanoseconds)
2.5
Figure 19: The top panels show the time evolution of Hy (left) and Ez (right) components at a selected point in the head. These solutions are obtained
M
using the fully tetrahedral DGTD-P2 method (solid line) and the hybrid DGTD-P2 Q2 (dashed line) and DGTD-P2 Q1 (dotted line) methods. In the bottom panels we plot the differences between the tetrahedral and hybrid DG solutions denoted as P2 − P2 Q2 and P2 − P2 Q1 .
464 465 466
d
463
References
[1] E. Abenius, U. Andersson, F. Edelvik, L. Eriksson, G. Ledfelt, Hybrid time domain solvers for the maxwell equations in 2d, Int. J. Numer. Meth. Eng. 53 (9) (2002) 2185–2199.
Ac ce pt e
462
[2] T. Apel, N. Düvelmeyer, Transformation of hexaedral finite element meshes into tetrahedral meshes according to quality criteria, Computing 71 (4) (2003) 293–304.
467
[3] C. A. Balanis, Advanced Engineering Electromagnetics, John Wiley & Sons, New York, 1989.
468
[4] K. Barrett, Multilinear jacobians for isoparametric planar elements, Finite Elem. Anal. Des. 40 (8) (2004) 821–853.
469
[5] M. Bergot, M. Duruflé, Higher-order discontinuous Galerkin method for pyramidal elements using orthogonal bases, Numer. Meth. Part. Diff.
470 471 472
Eq. 29 (1) (2013) 144–169.
[6] L. Branets, G. Carey, Extension of a mesh quality metric for elements with a curved boundary edge or surface, J. Comput. Inf. Sci. Eng. 5 (4) (2005) 302–308.
473
[7] K. Busch, M. König, J. Niegemann, Discontinuous Galerkin methods in nanophotonics, Laser & Photon. Rev. 5 (6) (2011) 773–809.
474
[8] G. Carey, Hexing the Tet, Communications in Numerical Methods in Engineering 18 (3) (2002) 223–227.
475
[9] L. Chew, Guaranteed-quality mesh generation for curved surfaces, in: 9th Annual ACM Symposium Computational Geometry, ACM Press,
476 477 478 479 480 481
1993, pp. 274–280.
[10] B. Cockburn, G. Karniadakis, C. Shu (eds.), Discontinuous Galerkin Methods. Theory, Computation and Applications, vol. 11 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin, 2000. [11] G. Cohen, X. Ferrieres, S. Pernet, A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell’s equations in time domain, J. Comput. Phys. 217 (2) (2006) 340–363. [12] R. Cools, Monomial cubature rules since “Stroud”: a compilation–part 2, J. Comput. Appl. Math. 112 (1-2) (1999) 21–27.
482
[13] R. Cools, An encyclopaedia of cubature formulas, J. Complex. 19 (3) (2003) 445–453, online database (http://nines.cs.kuleuven.be/ecf/).
483
[14] R. Cools, P. Rabinowitz, Monomial cubature rules since “Stroud”: a compilation, J. Comput. Appl. Math. 48 (3) (1993) 309–326.
26
Page 26 of 28
484 485 486 487
[15] S. Dosopoulos, B. Zhao, J.-F. Lee, Non-conformal and parallel discontinuous galerkin time domain method for maxwell’s equations: EM analysis of IC packages, J. Comput. Phys. 238 (0) (2013) 48–70. [16] D. Dunavant, Economical symmetrical quadrature rules for complete polynomials over a square domain, Int. J. Numer. Meth. Eng. 21 (10) (1985) 1777–1784.
488
[17] D. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle, Int. J. Numer. Meth. Eng. 21 (6) (1985) 1129–1148.
489
[18] D. Dunavant, Efficient symmetrical cubature rules for complete polynomials of high degree over the unit cube, Int. J. Numer. Meth. Eng.
491 492
23 (3) (1986) 397–407.
ip t
490
[19] C. Durochat, S. Lanteri, DGTD method on hybrid meshes for time domain electromagnetics, in: URSI International Symposium on Electromagnetic Theory (EMTS), 2010, pp. 992–995.
[20] C. Durochat, S. Lanteri, C. Scheid, High order non-conforming multi-element discontinuous Galerkin method for time-domain electromagnetics, in: International Conference on Electromagnetics in Advanced Applications (ICEAA), Cape Town, South Africa, 2012, pp. 379–382.
cr
493 494
[21] D. Eberly, 3D Game Engine Design, 2nd ed., Morgan Kaufmann Publishers, 2007.
496
[22] D. Eberly, Game Physics, 2nd ed., Morgan Kaufmann Publishers, 2010.
497
[23] F. Edelvik, G. Ledfelt, Explicit hybrid time domain solver for the Maxwell equations in 3D, J. Sci. Comput. 15 (1) (2000) 61–78.
498
[24] H. Fahs, Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic
503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521
an
502
[25] H. Fahs, High-order leap-frog based discontinuous Galerkin method for the time-domain Maxwell equations on non-conforming simplicial meshes, Numer. Math. Theor. Meth. Appl. 2 (3) (2009) 275–300.
[26] H. Fahs, Improving accuracy of high-order discontinuous Galerkin method for time-domain electromagnetics on curvilinear domains, Int. J. Comput. Math. 88 (10) (2011) 2124–2153.
M
501
[27] H. Fahs, Investigation on polynomial integrators for time-domain electromagnetics using a high-order discontinuous Galerkin method, Appl. Math. Model. 36 (11) (2012) 5466–5481.
[28] H. Fahs, L. Fezoui, S. Lanteri, F. Rapetti, Preliminary investigation of a non-conforming discontinuous Galerkin method for solving the
d
500
wave propagation, Int. J. Numer. Anal. Model. 6 (2) (2009) 193–216.
time-domain Maxwell equations, IEEE Trans. Magn. 44 (6) (2008) 1254–1257. [29] H. Fahs, A. Hadjem, S. Lanteri, J. Wiart, M.-F. Wong, Calculation of the SAR induced in head tissues using a high order DGTD method and
Ac ce pt e
499
us
495
triangulated geometrical models, IEEE Trans. Antennas Propag. 59 (12) (2011) 4669–4678. [30] H. Fahs, S. Lanteri, A high-order non-conforming discontinuous Galerkin method for time-domain electromagnetics, J. Comput. Appl. Math. 234 (4) (2010) 1088–1096.
[31] X. Ferrieres, J.-P. Parmantier, S. Bertuol, A. R. Ruddle, Application of a hybrid finite difference/finite volume method to solve an automotive EMC problem, IEEE Trans. Electromagn. Compat. 46 (4) (2004) 624–634. [32] G. Gassner, F. Lörcher, C.-D. Munz, J. Hesthaven, Polymorphic nodal elements and their application in discontinuous Galerkin methods, J. Comput. Phys. 228 (5) (2009) 1573–1590.
[33] P.-L. George, F. Hecht, E. Saltel, Automatic mesh generator with specified boundary, Comput. Methods Appl. Mech. Engrg. 92 (3) (1991) 269–288.
[34] S. Gottschalk, M. Lin, D. Manocha, OBB-tree: a hierarchical structure for rapid interference detection, in: Proceedings of ACM Siggraph’96, ACM, New York, USA, 1996, pp. 171–180.
[35] J. Hesthaven, From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex, SIAM J. Numer. Anal. 35 (2) (1998) 655–676.
522
[36] J. Hesthaven, C. Teng, Stable spectral methods on tetrahedral elements, SIAM J. Sci. Comput. 21 (6) (1999) 2352–2380.
523
[37] J. Hesthaven, T. Warburton, Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell’s equations, J. Comput.
524 525 526
Phys. 181 (1) (2002) 186–221. [38] J. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods: Algorithms, Analysis and Applications, Springer Texts Appl. Math., Springer-Verlag, Berlin, 2008.
27
Page 27 of 28
527
[39] P. Hubbard, Constructive solid geometry for triangulated polyhedra, Tech. Rep. CS-90-07, Brown University (1990).
528
[40] P. Knabner, S. Korotov, G. Summ, Conditions for the invertibility of the isoparametric mapping for hexahedral finite elements, Finite Elem.
529
Anal. Des. 40 (2) (2003) 159–172.
530
[41] P. M. Knupp, On the invertibility of the isoparametric map, Comput. Methods Appl. Mech. Engrg. 78 (3) (1990) 313–329.
531
[42] T. C. L. Zhang, H. Liu, A set of symmetric quadrature rules on triangles and tetrahedra, J. Comp. Math. 27 (1) (2009) 89–96.
532
[43] S. Mousavi, H. H. Xiao, N. Sukumar, Generalized gaussian quadrature rules on arbitrary polygons, Int. J. Numer. Meth. Eng. 82 (1) (2010)
536 537 538 539
ip t
535
[44] S. Pernet, X. Ferrieres, G. Cohen, High spatial order finite element method to solve Maxwell’s equations in time domain, IEEE Trans. Antennas Propag. 53 (9) (2005) 2889–2899.
[45] S. Piperno, M. Remaki, L. Fezoui, A nondiffusive finite volume scheme for the three-dimensional Maxwell’s equations on unstructured
cr
534
99–113.
meshes, SIAM J. Numer. Anal. 39 (6) (2002) 2089–2108.
[46] P. Ratiu, B. Hillen, J. Glaser, D. Jenkins, Medicine Meets Virtual Reality 11 - NextMed: Health Horizon, vol. 11, chap. Visible Human 2.0 the next generation, IOS Press, 2003, pp. 275–281.
us
533
540
[47] M. Rivero, F. Feito, Boolean operations on general planar polygons, Computers & Graphics 24 (6) (2000) 881–896.
541
[48] T. Rylander, A. Bondeson, Stability of explicit-implicit hybrid time-stepping schemes for Maxwell’s equations, J. Comput. Phys. 179 (2) (2002) 426–438.
an
542 543
[49] P. Schneider, D. Eberly, Geometric Tools for Computer Graphics, 1st ed., Morgan Kaufmann Publishers, 2003.
544
[50] S. Schnepp, T. Weiland, Efficient large scale electromagnetic simulations using dynamically adapted meshes with the discontinuous Galerkin
547 548
method, J. Comput. Appl. Math. 236 (18) (2012) 4909 – 4924.
[51] R. Sevilla, O. Hassan, K. Morgan, The use of hybrid meshes to improve the efficiency of a discontinuous galerkin method for the solution of maxwell’s equations, Comput. Struct.In Press.
M
545 546
[52] H. Songoro, M. Vogel, Z. Cendes, Keeping time with Maxwell’s equations, IEEE Microw. Mag. 11 (2) (2010) 42–49. [53] N. Sukumar, A. Tabarraei, Conforming polygonal finite elements, Int. J. Numer. Meth. Eng. 61 (12) (2004) 2045–2066.
550
[54] A. Taflove, S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed., Artech House, Norwood, MA, 2005.
Ac ce pt e
551
d
549
552
[55] S. Wandzura, H. Xiao, Symmetric quadrature rules on a triangle, Comput. Math. Appl. 45 (12) (2003) 1829–1840.
553
[56] R. B. Wu, T. Itoh, Hybrid finite-difference time-domain modeling of curved surfaces using tetrahedral edge elements, IEEE Trans. Antennas
554
Propag. 45 (8) (1997) 1302–1309.
28
Page 28 of 28
S0378-4754(14)00087-1 http://dx.doi.org/doi:10.1016/j.matcom.2014.04.002 MATCOM 4058
To appear in:
Mathematics and Computers in Simulation
Received date: Revised date: Accepted date:
19-8-2013 19-3-2014 7-4-2014
Please cite this article as: Hassan Fahs, High-order discontinuous Galerkin method for time-domain electromagnetics on non-conforming hybrid meshes, Mathematics and Computers in Simulation (2014), http://dx.doi.org/10.1016/j.matcom.2014.04.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript
High-order discontinuous Galerkin method for time-domain electromagnetics on non-conforming hybrid meshes Hassan Fahs∗
ip t
GPEM Group, Department of Materials & Structures, Nantes University, IFSTTAR, Route de Bouaye, 44344 Bouguenais Cedex, France
cr
Abstract
We present a high-order discontinuous Galerkin (DG) method for solving the time-dependent Maxwell equations
us
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
an
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
M
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
d
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
Ac ce pt e
methods allows for a notable reduction in computing time without sacrificing accuracy. Keywords: Maxwell’s equations, discontinuous Galerkin method, hybrid meshes, non-conformal meshes
1
1. Introduction
2
Nowadays, a variety of modeling strategies exist for the computer simulation of electromagnetic (EM) wave
3
propagations in the time-domain. In particular, the ability to deal with complex geometries has been extensively
4
studied in the context of unstructured methods as alternatives to the well established FDTD method [54]. In spite
5
of its versatility, the conventional FDTD method suffers from serious accuracy degradation when modeling curved
6
objects and small geometrical features. This is due to the commonly used staircase approximation (Cartesian grid)
7
which introduces large discretization errors, even when the grid size is very small. In contrast, unstructured solvers
8
such as the finite element time domain (FETD) [44] or the finite volume time domain (FVTD) [45] methods can
9
easily handle complex boundaries by dealing with unstructured tetrahedral meshes. Both methods are more expensive ∗ Corresponding
author Email address: [email protected], [email protected] (Hassan Fahs)
Preprint submitted to Elsevier
March 19, 2014
Page 1 of 28
than the FDTD method in terms of computing time and memory requirement. This motivates the development of
11
hybrid solution strategies, in which the computational efficiency of the Cartesian grid is combined with the geometric
12
flexibility of the unstructured mesh. The FDTD is then used in large homogeneous volumes of the computational
13
domain while the FVTD or FETD methods are employed in small regions near small scale features. For instance,
14
hybrid FETD-FDTD methods have been developed [56, 1], but they suffer from spurious reflections and late-time
15
instabilities, which limit their usefulness. A stable hybrid FETD-FDTD method has been described in [48]. To
16
achieve the time-stability, the scheme requires an implicit time integration in the unstructured elements but also a
17
transition layer of pyramids between the FDTD cubes and tetrahedra to ensure the mesh conformity. Low order
18
hybrid FDTD-FVTD methods are proposed in [23, 31]. These methods require the use of an intermediate layer of
19
hexahedra between the Cartesian grid and the tetrahedral mesh. Each cell in the layer is treated by the FVTD method,
20
which offers a correspondence between a quadrilateral face and two triangles in terms of the exchange of fluxes.
21
Again, weak instabilities have been reported for long-time simulations.
us
cr
ip t
10
The last few years have seen an increasing interest in so-called discontinuous Galerkin time domain (DGTD)
23
methods. The DG methods are a class of FE methods based on completely discontinuous piecewise polynomial
24
spaces for the numerical solution and the test functions. To obtain highly accurate and stable DG methods, suitable
25
numerical fluxes need to be designed over elemental interfaces. We refer to the lecture notes [10] and the textbook
26
[38] for details and history about DG methods. In particular, DG methods can easily handle elements of various
27
types and shapes [32], irregular non-conforming meshes [28, 15], and even locally varying polynomial degree [24],
28
and hence offer great flexibility in the mesh design. They also lead to (block-) diagonal mass matrices and therefore
29
yield fully explicit, inherently parallel methods when coupled with explicit time schemes. Furthermore, DG methods
30
are flexible with regards to the choice of the time-stepping scheme [27]. The DG methods have now acquired a
31
level of maturity and have successfully penetrated several scientific and technological communities following to their
32
adaptation to increasingly complex modeling contexts [29, 7, 50]. It is worthwhile to note that the DG method has
33
also been adopted for the first time in a commercial software as the time domain alternative of a very well known
34
EM wave simulation tool [52]. While DG methods on tetrahedral meshes have been extensively studied [37, 29, 15],
35
their development on hexahedral meshes has received relatively less attention [11]. One reason is that tetrahedral
36
meshing is highly automated while high quality hexahedral meshing of complex geometries commonly requires user
37
intervention and is labor intensive. This motivates the development of DG methods on conforming hybrid meshes
38
consisting of tetrahedra and hexahedra with transition layers of prisms or pyramids [5, 51]. Despite these advances,
39
the exploitation of non-conforming hybrid meshes in the context of a DG-based numerical methodology applied to
40
wave propagations, has been the subject of a very few studies so far. A 2D contribution in this context is the work
41
of Durochat et al. [19, 20] who proposed a DG method formulated on hybrid rectangular-triangular meshes for
42
Maxwell’s equations.
43
Ac ce pt e
d
M
an
22
We are concerned here with the design of a high-order nodal DGTD method for solving Maxwell’s equations 2
Page 2 of 28
on non-conforming hybrid meshes in 2D and 3D. The proposed method is based on centered numerical fluxes and
45
a fourth-order leap-frog time integration scheme. The hybrid mesh adopted here combines tetrahedral elements for
46
the discretization of irregularly shaped objects with hexahedral elements to fill the rest of the domain. Thus, the
47
resulting mesh is non-conformal in the sense that hanging nodes are allowed and no additional element types (e.g.,
48
pyramids or prisms) are needed to join tetrahedra with hexahedra. The non-conformal connection between these
49
elements is ensured thanks to the correct computation of flux terms by using efficient cubature rules over polygons.
50
The motivation behind this work stems partly from the following observations: first, the final mesh of an object is made
51
of separately constructed meshes and a conforming assembling is a very hard task. A good example to illustrate this
52
point is the interior penalty DGTD method formulated on non-conforming tetrahedral meshes presented in [15] for the
53
EM simulation of 3D integrated circuit packaging. For the complicated geometries involved in such applications, the
54
construction of a conforming mesh is a very hard task. Instead, a subdomain-based meshing approach is particularly
55
attractive and allowing meshes with different element-types is highly desirable in this context. Second, regarding the
56
computational effort, it turns out to be more efficient to compute on hexahedral meshes instead of tetrahedral meshes
57
to reach a desired error level. Therefore, the performance of the DG method can be increased by combining tetrahedra
58
with hexahedra. The outline of the paper is as follows. In Section 2 we briefly recall the basic features of DGTD
59
discretizations based on totally centered numerical fluxes and formulated on single element-type meshes composed
60
of triangles or quadrilaterals in 2D and tetrahedra or hexahedra in 3D. Then, we formulate the DG method on non-
61
conforming hybrid triangular-quadrilateral or tetrahedral-hexahedral meshes in Section 3. Subsequently, we analyze
62
the stability of the resulting algorithm and establish a sufficient CFL-like condition. Next, in Section 4, we discuss
63
some practical implementation issues including the generation of non-conforming hybrid meshes and the computation
64
of the flux matrices at hybrid interfaces between adjacent elements of different types. Finally, we present numerical
65
results for several test-cases in Section 5, and draw some conclusions of the study in Section 6.
66
2. Governing equations and DGTD formulation
68
cr
us
an
M
d
Ac ce pt e
67
ip t
44
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)
69
where ε and µ are respectively the permittivity and the permeability of the medium. These equations are posed on
70
a bounded domain Ω of Rd , d ≥ 2. At the domain boundary ∂Ω, we apply either perfect electric conductor (PEC)
71
conditions (n × E = 0) or absorbing boundary conditions (n × E + cµn × (n × H) = 0). Here n denotes the unit outward
72
normal to ∂Ω and c = (εµ)−1/2 is the speed of propagation.
73
We assume we dispose of a regular or irregular partition Th of Ω into ne non-overlapping elements κi , which can
74
be triangles or quadrilaterals in 2D and tetrahedra or hexahedra in 3D. Thus a (d −1)-dimensional face of each element 3
Page 3 of 28
75
76
77
78
κi in Th is allowed to contain hanging nodes. We shall suppose that each κi ∈ Th 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 ) P for all κi ∈ Th , where κr is either the open unit simplex T d = {ξξ = (ξ1 , . . . , ξd ) ∈ Rd : ξi ≥ 0, di=1 ξi ≤ 1} or the open unit hypercube Dd = {ξξ = (ξ1 , . . . , ξd ) ∈ Rd : 0 ≤ ξi ≤ 1} in Rd . 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 ) ∈ Nd0 , we define
80
spaces of polynomials of total degree ℓ ≥ 1 as follows
ip t
79
Qℓ (Dd ) = span{ξξα : 0 ≤ αi ≤ ℓ, 1 ≤ i ≤ d},
The dimension of each of these spaces depends on the spatial dimension d and on the degree ℓ and is given by d if κr = Dd , (ℓ + 1) N(ℓ, d) = (ℓ + d)!/ℓ!d! if κr = T d .
an
us
81
cr
α| ≤ ℓ}. Pℓ (T d ) = span{ξξα : 0 ≤ |α
Nℓ As a particular basis for Qℓ (Dd ) and Pℓ (T d ) consider the fundamental polynomials {L(ℓ) k }k=1 of multivariate Lagrange
83
interpolation, corresponding to a set of nodal points Sℓ = {η j , j = 1, . . . , Nℓ }, such that L(ℓ) k (η j ) = δ jk for η j ∈ Sℓ .
84
The distribution of the nodes η j is a key issue for the properties of the interpolation, especially for very high-order
85
approximations. For the unit simplex T d , we choose the nodal set derived from the electrostatics principle [35, 36],
86
while the Gauss-Legendre-Lobatto points are used for the hypercube Dd .
d
88
To each κi ∈ Th we assign an integer ℓi ≥ 1 that is the local polynomial degree inside κi ; collecting the ℓi and Ψκi in the vectors ℓ = {ℓi : κi ∈ Th } and Ψ = {Ψκi : κi ∈ Th }, respectively, we introduce the finite element space
Ac ce pt e
87
M
82
n d Vp (Ω, Th , Ψ ) = u ∈ [L2 (Ω)]d : u|κi ◦ Ψκi ∈ [Qℓi (κr )]d if Ψ−1 κi (κi ) = D ,
o d and u|κi ◦ Ψκi ∈ [Pℓi (κr )]d if Ψ−1 κi (κi ) = T , ∀κi ∈ Th ,
89
where [L2 (Ω)]d is the space of square-integrable functions on Ω. For two distinct elements κi and κk in Th , the
90
intersection aik = κi ∩ κk is a (d − 1)-dimensional face in Th , with unitary normal vector nik , oriented from κi towards
91
κk . For the boundary faces, the index k corresponds to a fictitious element outside the domain. Finally, we denote by
92
F I and F B the union of all interior and boundary faces of Th , respectively.
93
The derivation of the DG scheme for the simulation of electromagnetic wave propagation on purely simplicial
94
meshes has already been published in previous works [24, 30]. Therefore, we avoid the repetition of a comprehensive
95
derivation and refer the reader to these works for details. In brief, the electric and magnetic fields are approximated
96
inside each finite element κi by a linear combination of basis functions ϕi j = L(ℓj i ) ◦ Ψ−1 κi of degree ℓi with support κi
97
and with time-dependent coefficient functions Ei j (t) and Hi j (t) as follows Ei (x, t) =
X
Ei j (t)ϕi j (x),
Hi (x, t) =
X
Hi j (t)ϕi j (x).
(2)
1≤ j≤Ni
1≤ j≤Ni
4
Page 4 of 28
Here, the index j indicates the jth basis function and Ni = N(ℓi , d) is the number of degrees of freedom inside κi .
99
¯ i , respectively, the column vectors (Ei j )1≤ j≤Ni and (Hi j )1≤ j≤Ni . As usual for DG schemes, the We denote by E¯ i and H
100
Maxwell system, Eq. (1), is multiplied by a test function ϕ ∈ span{ϕi j , 1 ≤ j ≤ Ni } and integrated over each κi .
101
After integration by parts, inserting the DG approximation, Eq. (2), and after applying a centered numerical flux in
102
the boundary integrals, the semi-discrete formulation of the scheme in the physical element κi reads as =
¯i− Ki H
X
¯k− Sik H
X
¯k + Sik E
aik ∈F I
¯i Mµi ∂t H
= −Ki E¯ i +
¯ k, Sik H
X
Sik E¯ k ,
aik ∈F B
aik ∈F I
aik ∈F B
(3)
where the mass matrices Mγi (γ stands for ε or µ), the stiffness matrix Ki and the interface (or flux) matrix Sik write Z Z Z 1 1 ϕi j · ϕil , (Ki ) jl = (Mγi ) jl = γi ϕi j · ∇ × ϕil + ϕil · ∇ × ϕi j , (Sik ) jl = ϕi j · (ϕkl × nik ). (4) 2 κi 2 aik κi
us
103
X
cr
¯i Mεi ∂t E
ip t
98
These matrices involve volume and surface integrals which are evaluated numerically over the master element, κr ,
105
using the mapping Ψκi . For straight-sided simplices and parallelepipedic elements, all integrals can be computed using
106
affine mappings. Since the Jacobian of the affine mapping is constant, the matrices in Eq. (4) can be precomputed and
107
stored for κr in advance of the main calculation once and for all. Straight-sided convex quadrilateral and hexahedral
108
elements are expressed by isoparametric bilinear and trilinear mappings [41, 40]. These mappings are not linear and
109
the Jacobian is not constant over a quadrilateral and hexahedron anymore. Thus, all metric terms relating to Ψκi are
110
evaluated, and need to be stored for each node in the elements. Finally, we also allow for curved elements which are
111
usually present inside Th (e.g., around a boundary layer or a curved material interface) or in the vicinity of ∂Ω. These
112
are introduced in order to avoid geometrical errors and to achieve optimal convergence for high-order polynomials.
113
In our implementation we use quadratic interpolation of curved boundaries, and define isoparametric mappings based
114
on that approximation of the surface [26, 4, 6]. For such mappings, the Jacobians are not constants and all matrices
115
should be computed and stored for each curved element.
Ac ce pt e
d
M
an
104
116
Once the local mapping Ψκi is introduced for each κi , the integrals in Eq. (4) are evaluated using cubature formulas
117
which are able to accurately integrate high-order polynomials over the master element κr . A cubature is a set of N c
118
k=Nc k=Nc with N c associated weights {ωck }k=1 , where the number of points N c depends on the desired 2D or 3D points {λck }k=1
119
maximum order of polynomial required to be accurately integrated. Surveys of cubature formulas were performed
120
in [14, 12] and an on-line database containing many of the best known rules is described in [13]. Here, we adopt
121
the efficient high-order symmetric rules proposed in [42, 55, 17, 16, 18] for hypercubes and simplices. Finally, edge
5
Page 5 of 28
integrals are calculated using Gauss-Legendre quadrature. We evaluate γ
(Mi ) jl
=
γi
k=N Xc k=1
(Ki ) jl
=
(Sik ) jl
=
h i i) , ωck L(ℓj i ) L(ℓ kJ k κ i l c λk
k=N i 1 Xc c h (ℓi ) ωk L j (jκi ) jl ∇ × Ll(ℓi ) + Ll(ℓi ) (jκi ) jl ∇ × L(ℓj i ) c , λk 2 k=1 k=N i 1 Xc c h (ℓi ) (ℓi ) ωk L j Ll k∇Jκi k c , λk 2 k=1
(5)
ip t
122
where each term in the brackets of the right hand sides are evaluated at the cubature points, kJκi k is the determinant of
124
the Jacobian matrix, and jκi = kJκi kJ−1 κi .
127
us
126
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 Pe (resp. H) of size Ng = ni=1 Ni . Then system (3) can be rewritten as Mε ∂t E
=
Mµ ∂t H
= −KE + AE − BE,
KH − AH − BH,
an
125
cr
123
(6)
where we have the following definitions and properties: (a) Mε , Mµ and K are Ng × Ng block diagonal matrices with
129
diagonal blocks equal to Mεi , Mi , and Ki respectively. Mε and Mµ are symmetric positive definite matrices, and K is
130
a symmetric matrix; (b) A is also a Ng × Ng symmetric block sparse matrix, whose non-zero blocks are equal to Sik
131
when aik ∈ F I ; (c) B is a Ng × Ng skew-symmetric block diagonal matrix, whose non-zero blocks are equal to Sik
132
when aik ∈ F B . Let S = K − A − B; the system (6) rewrites as
Ac ce pt e
d
µ
M
128
Mε ∂t E
=
SH,
Mµ ∂t H
=
−S† E.
(7)
133
where the symbol † stands for the matrix transpose. To define a fully discrete scheme, we divide the time interval
134
[0, T ] into Mt uniform subintervals such that 0 = t0 < t1 < · · · < t Mt = T , where tn = n∆t and ∆t is the fixed time step.
135
Time integration of the semi-discrete system, Eq. (7), is done using a fourth-order explicit leap-frog scheme. Let us
136
137
denote En = E(tn ) and Hn+1/2 = H(tn+1/2 ), the fully discrete central DGTD method writes En+1 − En 1 Mε = GHn+ 2 , ∆t (8) 1 n+ 23 − Hn+ 2 † n+1 Mµ H = −G E , ∆t 2 µ −1 † ε −1 where G = S I − ∆t 24 (M ) S (M ) S and I is the identity matrix. The resulting DGTD method is analyzed in
138
[30, 25] on purely simplicial meshes where it is shown that the method is non-dissipative, conserves a discrete form
139
of the EM energy and is stable under the CFL-like condition 1
1
∆t ≤ 2/k(Mε )− 2 G(Mµ )− 2 k, 1
140
where k · k denotes any matrix norm, and the matrix (Mγ )− 2 is the inverse square root of Mγ (γ stands for ε or µ). 6
Page 6 of 28
141
3. DGTD method on hybrid tetrahedral-hexahedral meshes
142
3.1. The numerical scheme The set of elements κi of the mesh is now assumed to be partitioned into two subsets: one made of the tetrahe-
144
dral (triangular in 2D) elements and the other one gathering the hexahedral (quadrilateral in 2D) elements. In the
145
following, these two subsets are respectively referred as Sτ and Sh . There is no need of a particular assumption on
146
the connectivity of the two subsets. In order to further describe this scheme, we introduce additional definitions. First,
147
the problem unknowns are reordered as E τ E = E
us
h
where subvectors with a τ subscript (resp., a h subscript) are associated to the elements of the set Sτ (resp., the set Sh ). We deduce from this partitioning of the unknown vectors the following decompositions of the system matrices B K Mµ 0 Mε 0 0 0 τ τ τ τ . , B = , K = , Mµ = Mε = 0 Mε 0 B 0 K 0 Mµ h
h
151
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 Sik is decomposed as A ττ Aτh A = , A A
154
155
hh
where Aττ and Ahh are symmetric matrices, and Aτh = A†hτ . Introducing the two matrices Sτ = Kτ − Aττ − Bτ and
Ac ce pt e
153
d
hτ
152
h
h
M
150
an
149
cr
H τ and H = , H
h
148
ip t
143
Sh = Kh − Ahh − Bh , the global system of ordinary differential equations (7) can be split into two systems ε Mτ ∂t Eτ = Sτ Hτ − Aτh Hh , Mµτ ∂t Hτ = −S†τ Eτ + Aτh Eh , ε Mh ∂t Eh = Sh Hh − Ahτ Hτ , Mµh ∂t Hh = −S†h Eh + Ahτ Eτ . Then, the proposed DGTD scheme on hybrid meshes can be written as ! n+1 − Enτ n+ 1 n+ 1 ε Eτ M = Gτ Hτ 2 − Gτh Hh 2 , τ ∆t n+ 3 1 H 2 − Hn+ 2 µ τ τ = −G†τ En+1 + Gτh En+1 τ h , Mτ ∆t n+1 Eh − Enh ε = Mh ∆t n+ 3 1 H 2 − Hn+ 2 µ h h = M h ∆t
n+ 12
Gh Hh
n+ 21
− Ghτ Hτ
(9)
(10a)
, (10b)
−G†h En+1 h
+
Ghτ En+1 τ ,
7
Page 7 of 28
156
where
! ∆t2 µ −1 † ε −1 Gτ = Sτ I − (Mτ ) Sτ (Mτ ) Sτ , 24 ! ∆t2 µ −1 † ε −1 (Mh ) Sh (Mh ) Sh , Gh = Sh I − 24
! ∆t2 µ −1 † ε −1 I− (Mh ) Aτh (Mτ ) Aτh , 24
Gτh = Aτh Ghτ = G†τh .
The DGTD method on triangular-quadrilateral meshes or tetrahedral-hexahedral meshes, Eq. (10), is referred to as
158
DGTD-P p Qq in the following. If the mesh contains fully triangular or tetrahedral elements (resp. quadrilateral or
159
hexahedral elements), then Eq. (10a) (resp. Eq. (10b)) is referred to as DGTD-P p method (resp. DGTD-Qq method).
160
3.2. Energy conservation
cr
We define the following discrete electromagnetic energy En =
us
161
ip t
157
1 1 1 n † ε n n− 1 n− 1 µ n+ µ n+ (Eτ ) Mτ Eτ + (Hτ 2 )† Mτ Hτ 2 + (Enh )† Mεh Enh + (Hh 2 )† Mh Hh 2 . 2
(11)
In the following, we shall verify that En is exactly conserved in absence of absorbing boundary conditions, i.e.,
163
∆E := En+1 − En = 0.
164
Lemma 1. Using the hybrid DGTD-P p Qq method, Eq. (10), for solving the Maxwell system, Eq. (1), with PEC
165
boundary conditions only, the energy En is exactly conserved.
166
n Proof. We denote by Eτ/h2 = (Enτ/h + En+1 τ/h )/2. The variation of the energy E through one time step is given by
M
an
162
=
1 n+ 1 n+ 3 n− 1 µ ) Mετ [En+1 − Enτ ] + (Hτ 2 )† Mτ [Hτ 2 − Hτ 2 ] + τ 2 1 n+ 1 n+ 3 n− 1 n+ 21 † ε n+1 n (Eh ) Mh [Eh − Eh ] + (Hh 2 )† Mµh [Hh 2 − Hh 2 ] 2 n+ 12 †
(Eτ
Ac ce pt e
∆E
d
n+ 1
=
n+ 21 †
n+ 12
n+ 21
n+ 21
∆t(Eτ
) Gτ Hτ
∆t(Eh )† Gh Hh
=
n+ 21
n+ 1
n+ 12
) [Gτ − Gτ ]Hτ
∆t(Eh 2 )† [Gh − Gh ]Hh
167
which completes the proof.
168
3.3. Stability analysis
169
170
n+ 12
n+ 12
n+ 12
) Gτh Hh
− ∆t(Eh )† Ghτ Hτ
n+ 21 †
∆t(Eτ
n+ 21 †
− ∆t(Eτ
n+ 12 †
n+ 12
n+ 21
n+ 12
− ∆t(Hτ
) G†τ Eτ
− ∆t(Hh )† G†h Eh
n+ 12 †
n+ 21
n+ 1
n+ 1
− ∆t(Eτ
) [Gτh − G†hτ ]Hh
n+ 12 †
n+ 21
n+ 21
n+ 12
+ ∆t(Hτ
) Gτh Eh
+
+ ∆t(Hh )† Ghτ Eτ
+
− ∆t(Eh 2 )† [Ghτ − G†τh ]Hτ 2 ,
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.
8
Page 8 of 28
171
Lemma 2. For the hybrid DGTD-P p Qk method, Eq. (10), the electromagnetic energy, En , is a positive definite
172
quadratic form of all numerical unknowns provided that
1 1 , βτ + βτh βh + βτh
!
with
h
1
Proof. Using the second relations of Eq. (10a) and Eq. (10b), the energy En , Eq. (11) can be written as 1 n † ε n n− 1 n− 1 n− 1 n− 1 (Eτ ) Mτ Eτ + (Hτ 2 )† Mµτ Hτ 2 − ∆t(Hτ 2 )† G†τ Enτ − ∆t(Hτ 2 )† Gτh Enh En = 2 n− 1 n− 1 n− 1 n− 1 +(Enh )† Mεh Enh + (Hh 2 )† Mµh Hh 2 − ∆t(Hh 2 )† G†h Enh − ∆t(Hh 2 )† Ghτ Enτ .
us
174
where k.k denotes any matrix norm, and the matrix (Mγτ/h )− 2 is the inverse square root of Mγτ/h (γ stands for ε or µ).
cr
173
(12)
ip t
∆t ≤ 2 min
1 1 βτ = k(Mετ )− 2 Gτ (Mµτ )− 2 k, 1 µ 1 βh = k(Mεh )− 2 Gh (Mh )− 2 k, βτh = k(Mµτ )− 21 Gτh (Mε )− 21 k,
The mass matrices Mετ/h and Mµτ/h are symmetric positive definite and we can construct their square root (also sym-
176
metric positive definite) denoted by (Mετ/h ) 2 and (Mµτ/h ) 2 , respectively. We deduce from Eq. (11) that
1
≥
1
1 1 1 1 1 1 1 n− 1 n− 1 n− 1 k(Mετ ) 2 Enτ k2 + k(Mµτ ) 2 Hτ 2 k2 − ∆tβτ k(Mµτ ) 2 Hτ 2 k k(Mετ ) 2 Enτ k − ∆tβτh k(Mµτ ) 2 Hτ 2 k k(Mεh ) 2 Enh k 2 1 1 1 1 1 1 n− 1 n− 1 n− 1 +k(Mεh ) 2 Enh k2 + k(Mµh ) 2 Hh 2 k2 − ∆tβh k(Mµh ) 2 Hh 2 k k(Mεh ) 2 Enh k − ∆tβhτ k(Mµh ) 2 Hh 2 k k(Mετ ) 2 Enτ k ,
M
En
an
175
1
≥
d
En
178
1
where βτ , βh and βτh are defined in Eq. (12), and βhτ = k(Mµh )− 2 Ghτ (Mετ )− 2 k. Then estimate ab ≤ (a2 + b2 )/2 leads to 1 1 1 1 ∆t ∆t 1 1 1 µ 1 n− µ 1 n− µ 1 n− k(Mετ ) 2 Enτ k2 + k(Mτ ) 2 Hτ 2 k2 − βτ k(Mτ ) 2 Hτ 2 k2 + k(Mετ ) 2 Enτ k2 − βτh k(Mτ ) 2 Hτ 2 k2 + k(Mεh ) 2 Enh k2 2 2 2 1 1 ∆t ∆t 1 1 1 1 1 1 n− n− n− 1 +k(Mεh ) 2 Enh k2 + k(Mµh ) 2 Hh 2 k2 − βh k(Mµh ) 2 Hh 2 k2 + k(Mεh ) 2 Enh k2 − βhτ k(Mµh ) 2 Hh 2 k2 + k(Mετ ) 2 Enτ k2 . 2 2
Ac ce pt e
177
We then sum up the lower bounds for the En to obtain βτ + βhτ 1 1 1 − ∆t k(Mετ ) 2 Enτ k2 + En ≥ 2 2 βh + βτh 1 1 1 − ∆t k(Mεh ) 2 Enh k2 + 2 2
βτ + βτh 1 1 n− 1 1 − ∆t k(Mµτ ) 2 Hτ 2 k2 + 2 2 1 βh + βhτ 1 n− 1 1 − ∆t k(Mµh ) 2 Hh 2 k2 . 2 2
179
Since the parameters ε and µ are piecewise constant and Gτh = G†hτ , we have that βhτ = βτh . Then, under the condition
180
proposed in Lemma 2, the energy En is a positive definite quadratic form of the numerical unknowns Enτ/h and Hτ/h2 ,
181
and Eq. (12) states a sufficient condition for the stability of the hybrid DGTD-P p Qq method.
182
4. Implementation issues
183
n− 1
For generality, we shall limit the discussions to the 3D case and regard the 2D case as a natural simplification.
9
Page 9 of 28
184
4.1. Nonconforming hybrid mesh generation A simple way for generating a non-conforming hybrid tetrahedral-hexahedral mesh without pyramidal or prismatic
186
elements is based on element-type conversion. This indirect approach starts with the generation of a conforming or
187
non-conforming mesh with only one type of elements (pure tetrahedra or pure hexahedra). It is partially converted
188
to hexahedral or tetrahedral mesh by applying a conventional subdivision scheme. For instance, a tetrahedron can
189
be split into four hexahedra [8] while a hexahedron can be decomposed into five or six tetrahedra [2]. The resulting
190
hybrid mesh is non-conformal where the hybrid interface between a tetrahedron and a hexahedron is either a triangle
191
or a quadrilateral, see Figure 1. Although this approach avoids some difficulties, it rapidly increases the number of
192
elements and tends to introduce badly shaped elements which may deteriorate their numerical properties. We will not
193
pursue this approach further here, and will propose an alternative method.
us
M
Neighbor hexahedron
Hybrid interface
an
Hexahedron converted to 5 tetrahedra
Tet
ing
rm
.
fo on −c s on ace d n erf bri int
..
.
0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1111111111111 0000000000000 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 10 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 1010 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 00 11 10 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111 0000000000000 1111111111111
Hy
rah e to 4dron c hex onv ahe erte dra d
Neighbor tetrahedron
Hybrid interface
111111 000000 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 000000 111111 000000 111111 0000000 1111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 0000000 1111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111
cr
ip t
185
Figure 1: Examples of hybrid interfaces when the hybrid tetrahedral-hexahedral mesh is generated using the element-type conversion approach. As
d
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
194
195
Ac ce pt e
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 S domain Ω is decomposed into subdomains Ωi such that Ω = i Ωi . Then, an individual mesh for each subdomain Ωi
196
is generated using a suitable meshing algorithm and the most appropriate element type. The resulting meshes should
197
not share nodes, edges, faces, or cells. Finally, we combine these submeshes to obtain the non-conformal hybrid
198
mesh of the whole domain Ω. This strategy of non-conforming hybrid mesh generation has several advantages over
199
other approaches based on hybrid conforming mesh. For instance, it allows one (i) to avoid the intermediate layer
200
of pyramids as transition elements between triangular and quadrilateral faces; (ii) to relax the node-to-node matching
201
requirement at the interface between the meshes; (iii) to change or replace certain portions of the domain without
202
changing the entire mesh. These features are particularly useful for parametric modelling of complex geometries
203
where the object is built from several parts designed separately. Figure 2(a)-(b) illustrates examples of hybrid meshes
204
generated by this strategy.
205
After combining the meshes, we need to find the hybrid non-conforming interfaces between triangular and quadri-
206
lateral faces in order to compute the flux matrices, see Figure 2(c). First, we locate the boundary faces of each
207
subdomain mesh and store them in a list including their element and face location. Next, we use the fact that only 10
Page 10 of 28
the boundary faces of two subdomains that shared boundaries may participate in a partial face-face connection. Let
209
us assume that two subdomains Ωi and Ω j share a boundary. For each of the boundary faces of the mesh of Ωi , we
210
perform two tests: the first is an intersection test with every member of the boundary face list of the mesh of Ω j . This
211
just indicates whether or not an intersection exists. The second is a find-intersection test that involves computing the
212
set of intersections between triangular and quadrilateral faces. To perform these intersection tests, we use the class
213
Intersection provided by the Eberly’s Wild Magic library [22]. For the reader’s convenience, we briefly describe
214
the general idea of this class. The test-intersection is based on the method of separating axes [34], a method for deter-
215
mining whether or not two convex objects are intersecting. A test for nonintersection of two convex objects is simply
216
stated: if there exists a line for which the intervals of projection of the two objects onto that line do not intersect, then
217
the objects do not intersect. Such a line is called a separating axis. For a pair of convex polygons in 3D, only a finite
218
set of direction vectors needs to be considered for separation tests. The polygons normals are potential separating
219
directions. If either normal direction separates the polygons, then no intersection occurs. However, if neither normal
220
direction separates the polygons, then two possibilities exist for the remaining potential separating directions. The
221
first case is that the polygons are coplanar and the remaining potential separating directions are those vectors in the
222
common plane and perpendicular to the polygon edges. The second case is that the polygon planes are not parallel
223
and that each polygon is split by the plane of the other polygon. The remaining potential separating directions are
224
the cross products of two edges, one from each polygon. In summary, for the intersection between a triangular and
225
a quadrilateral face in 3D, a maximum of 19 separating directions need to be tested in the worst case. Given a line,
226
it is straightforward to project the vertices of a polygon onto the line and compute a bounding interval for those pro-
227
jections. The two intervals obtained from the two polygons are easily compared for overlap. If an intersection exist,
228
we perform a find-intersection query using Boolean operations on polygons [47, 39]. Here the intersection between
229
a triangular and a quadrilateral face (both lies in the same plane) is a convex polygon with at most seven vertices.
230
Let Q be a quadrilateral and T be a triangle with vertices (v1 , v2 , v3 ) that are counterclockwise ordered. We define
231
half-planes, πi j , 1 ≤ i < j ≤ 3, bounded by the lines (vi v j ) and containing the third vertex of the triangle T . The
232
intersection between T and Q is defined as T ∩ Q = ((Q ∩ π12 ) ∩ π23 ) ∩ π13 . The problem is reduced to find the
233
intersection between a polygon and a half-plane which is relatively straightforward. For a more complete exposition
234
of the theory of intersection tests, we refer the reader to the books [49, Chap. 7 & 11]-[21, Chap. 15].
235
4.2. Non-conforming hybrid interface matrix calculation
Ac ce pt e
d
M
an
us
cr
ip t
208
236
A key ingredient for the implementation of the DG method on hybrid meshes is the computation of the matrix
237
Aτh involving integrals over a hybrid interface. We concentrate here on the evaluation of such integrals. In fact, the
238
other matrices in Eq. (9) which involve integrals over an element or an interface between two elements of the same
239
type are evaluated as specified in Eq. (5). We just have to utilize the correct set of basis functions corresponding to
240
the particular element type.
241
Let aik = τi ∩ hk be a hybrid interface between a tetrahedron τi and a hexahedron hk . Here aik is a n-gon with 11
Page 11 of 28
1111111 0000000 ■ 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 0000000 1111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111
.
Hy
.
fo on −c on es d n ac bri interf
■
ip t
ng
(c)
cr
(b)
..
.
i rm
(a)
111111 000000 000000 111111 000000 111111 000000 111111 000000 111111 00000000000000 11111111111111 000000 111111 00000000000000 11111111111111 000000 111111 00000000000000 11111111111111 000000 111111 00000000000000 11111111111111 000000 111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 ■11111111111111 00000000000000 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 00000000000000 11111111111111 ■
Figure 2: Examples of non-conformal hybrid meshes generated using the combining meshes strategy. In (a) we show a triangular-quadrilateral
us
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
an
node (filled circle).
242
n = 3, 4, 5, 6, 7, that is, aik is a convex polygon with n edges (see Figure 3). Let also R be a regular polygon in R2 and
243
Υ : R ∋ ξ → x ∈ aik ⊂ R3 be an explicit parametrization of aik . The matrix Aτh is block diagonal, whose non-zero blocks are equal to the interface matrix Sik when aik is a hybrid interface. The matrix Sik is defined as Z Z ) −1 k i) ξ ) L(q ξ ) k∇Υ(ξξ )kdξξ L(p ◦ Ψ−1 θi j (x)ϑkl (x)dx = (Sik ) jl = hk ◦ Υ(ξ j ◦ Ψτi ◦ Υ(ξ l R
aik
M
244
(13)
where θi j ∈ P pi (τi ) and ϑkl ∈ Qqk (hk ) are the basis functions of degrees pi and qk in the elements τi and hk , respec-
246
tively. The integrals in Eq. (13) are calculated numerically and stored for each hybrid interface aik . To this end, two
247
approaches are possible
Ac ce pt e
d
245
248
1. Partitioning the interface aik (a n-gon) into n triangles and then performing cubature rules over each triangle
249
[55, 17]. For n-gons with n ≥ 4, we use the centroid of aik to partition it into triangles. The matrix Sik is
250
computed with the following formula
(Sik ) jl =
251
252
253
254
θi j (x)ϑkl (x)dx =
aik
s=n Z X s=1
θi j (x)ϑkl (x)dx,
(14)
∆s
where ∆ s is a triangle formed by two successive vertices of aik and its centroid as illustrated in Figure 4(a). The P parametrization of ∆ s can be defined as x = Υ(ξξ ) = 3i=1 φi (ξξ )vi , where ξ = (ξ1 , ξ2 ) ∈ T 2 , φi are the barycentric coordinates in T 2 and vi denotes the three vertices of ∆ s . To compute the last integral in Eq. (14), we lay down Nc cubature points λn ∈ T 2 and weights ωn . This is achieved by Z
∆s
255
Z
θi j (x)ϑkl (x)dx =
Nc X
(p )
(q )
L j i (λτn )Ll k (λhn )ωn k∇Υ(λn )k,
(15)
n=1
3 h −1 3 where λτn = Ψ−1 τi ◦ Υ(λn ) ∈ T and λn = Ψhk ◦ Υ(λn ) ∈ D .
12
Page 12 of 28
2. An alternative is to use cubature rules over polygons. However, the design and development of such rules has
257
not yet reached a mature stage. To our knowledge this approach has not been used previously in the context
258
of DG methods. Recently, Mousavi et al. [43] presented an algorithm to build symmetric cubature rules for
259
arbitrary regular polygons. This algorithm combines elements of group theory and numerical optimization and
260
was originally proposed by Wandzura and Xiao [55] for the construction of cubature rules over triangles. Let Pn
261
be a reference regular n-gon with vertices pm = (cos 2πm/n, sin 2πm/n) for m = 1, 2, . . . , n in the ξ -coordinate
264
point p ∈ Pn defined as [53]
αi (ξξ ) , φi (ξξ ) = Pn ξ) j=1 α j (ξ
α j (ξξ ) =
(1 − |ξξ |2 ) sin3 (π/n) cos(π/n) , Ai−1 Ai
cr
263
system, see Figure 4(b). The parametrization Υ from Pn to aik is obtained via the isoparametric mapping P x = Υ(ξξ ) = ni=1 φi (ξξ )vi , where vi are the vertices of the n-gon aik and φi are the Laplace shape functions at ξ = (ξ1 , ξ2 ) ∈ Pn ,
us
262
ip t
256
where αi (ξξ ) is the Laplace weight function, Ai is the area of the triangle [p, pi , pi+1 ] and Ai−1 is the area of the
266
triangle [p, pi−1 , pi ]. The Laplace shape functions are non-negative, interpolate nodal data and form a partition
267
of unity. They are also linear on the polygon boundaries which ensures that the parametrization Υ(ξξ ) is an affine
268
mapping. Once Υ(ξξ ) is given explicitly, the integral in Eq. (13) is evaluated by laying down Nc cubature points
269
on Pn and proceed as in Eq. (15).
M
an
265
We compare the two approaches by evaluating the integral of x6 y6 z6 over the regular pentagon, hexagon and heptagon.
271
The relative errors between exact and approximate values are shown in Figure 5 as a function of the number of cubature
272
points. We observe that for a relative error of 10−8 the polygonal cubatures require 3 to 5 times fewer points than the
273
partitioning approach. Hence, for the numerical results given in Section 5, we adopt the second approach since it is
274
more practical and requires the least number of cubature points for a given degree of accuracy.
Ac ce pt e
d
270
●
●
●
aik
aik
●
●
●
●
aik
●
●
●
● ● ●
● ●
●
●
aik
● ●
●
●
aik
● ●
● ●
Figure 3: Five possible configurations for the hybrid interface aik between a quadrilateral face and a triangular face in 3D.
275
5. Numerical results
276
5.1. Two-dimensional examples
277
We consider here the 2D TMz -polarized Maxwell equations, that is, we solve Eq. (1) for the fields H = (H x , Hy , 0)
278
and E = (0, 0, Ez). Two examples with analytical solutions are studied in detail to illustrate the accuracy and per-
279
formance of the DG method on hybrid triangular-quadrilateral meshes. The exact time-domain solutions of these
280
problems can be found in [24]. 13
Page 13 of 28
. . ξ2
p2
.
T2
ξ1
.
v3
x1
.
3
v2
∆s
.
■
.
x3
.
p
.
x2
1
.
p6
ξ1
.
. . . . P6
p4
.
v1
p
p
v2
v3
5
x2
x1
(a)
v1
aik
x3
v4
.
ip t
ξ2
.
.
v6
cr
. .
v5
(b)
us
Figure 4: Two numerical integration schemes used for the evaluation of integrals over a polygon aik . In (a) we show a scheme based on the partition of aik into triangles ∆ s . The scheme illustrated in (b) consists of using cubature rules over regular polygons. 100 Hexagonal cubature Partitioning approach
10-4 135 points
198 points
10-8 47 points 0
50
100 150 200 Number of cubature points
47 points
46 points 250 0
50
100 150 200 Number of cubature points
d
10-10
231 points
M
Relative error
10-2
10-6
Heptagonal cubature Partitioning approach
an
Pentagonal cubature Partitioning approach
250 0
50
100 150 200 Number of cubature points
250
Ac ce pt e
Figure 5: Convergence curves of two numerical integration approaches used to evaluate the integral of x6 y6 z6 over the regular pentagon (left), hexagon (center) and heptagon (right).
281
5.1.1. Eigenmode in a PEC cavity
282
In the first example, we propose to validate the implementation of the proposed DGTD-P p Qq method on non-
283
conforming hybrid meshes by performing a numerical convergence test. To this end, we consider the propagation of
284
an eigenmode which is a standing wave of frequency f = 212 MHz and wavelength of 1.4 m in a unitary vacuum-filled
285
square cavity with PEC walls. The computations are performed on three successively refined hybrid meshes, named
286
MSi for i = 1, 2, 3. Each mesh MSi consists of 5 × 22i−1 triangles and 3 × 22i−1 quadrilaterals, as shown in Figure 6.
287
We present results with interpolation orders, p and q, up to degree 8 such that p = q. For each interpolation order,
288
we use the maximum time step, ∆t, allowed by the stability condition of the hybrid method. In Figure 7 (left) we
289
illustrate the h-convergence graphs as a function of the total number of degrees of freedom (#DOF). As expected, for
290
smooth solutions the convergence rate behaves as O(∆t4 + h p ) even for non-conforming hybrid meshes. Also shown in
291
Figure 7 (right) is the p-convergence rate in linear-log scale, observing an exponential decay of the error as exp(−s p p),
292
with decay constant s p = log[E(p − 1, q − 1)/E(p, q)] ≃ 2, where E(p, q) is the L2 error when the interpolation orders
293
p and q are used. By inspecting the results in Figure 7, we observe that the gain in accuracy is notable when refining 14
Page 14 of 28
the mesh or increasing the interpolation orders. Furthermore, to achieve a certain accuracy level, the DG methods
295
require less #DOF for higher interpolation orders, as typically found with other DG methods based on fully simplicial
296
meshes [25, 30]. 1
0.75
0.75
0.75
0.5
0.5
0.5
0.25
0.25
0.25
0
0 0
0.25
0.5
0.75
cr
1
us
1
ip t
294
0
1
0
0.25
0.5
0.75
1
0
0.25
0.5
0.75
1
an
Figure 6: Sequences of three successively refined hybrid meshes used for assessing convergence of the eigenmode problem.
100
10-2
10-2
p=q=2, rate=2.08 p=q=3, rate=2.92
10-6
L2 error
10-4
10-4
p=q=4, rate=3.87
d
L2 error
M
p=q=1, rate=1.25
p=q=5, rate=4.10
10-8
10-6
Ac ce pt e
p=q=6, rate=4.24
p=q=7, rate=4.08
p=q=8, rate=4.13
10-10
10
10-8 2
100
(DOF)1/2
3
4
5 6 Interpolation order
7
8
Figure 7: On the left is shown the L2 -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(∆t4 + h p ) for all orders. On the right, we plot the slope of the error for the mesh MS1 , confirming the spectral convergence of the hybrid DGTD-P p Qq method.
297
5.1.2. Plane-wave scattering by a dielectric circular cylinder
298
As an illustration of the capability of the hybrid DG method to handle materials, let us consider plane wave
299
scattering by a penetrable circular cylinder with a radius of 0.6 m consisting of a dielectric nonmagnetic material
300
with a relative permittivity εr = 12.0. The frequency of the incident wave is f = 300 MHz and the wavelength is
301
λ = 0.288 m. The surrounding medium is assumed to be vacuum. The computational domain Ω is bounded by a
302
square of side length 3.2 m centered at (0, 0). A Silver-Müller absorbing condition is applied on the boundary of
303
the square. We compare the performance of the DGTD-P p Qq method using a non-conforming hybrid mesh with the
304
DGTD-P p method on a fully triangular mesh, see Figure 8. The hybrid mesh is generated by combining submeshes 15
Page 15 of 28
of three subdomains, Ω1 = [−0.35, 0.35]2 m2 (inside the cylinder), Ω2 = [−0.7, 0.7]2 m2 \Ω1 (close to cylinder and
306
completely enclosing it) and Ω3 = [−1.6, 1.6]2 m2 \(Ω1 ∪ Ω2 ) such that Ω = Ω1 ∪ Ω2 ∪ Ω3 . The regions Ω1 and
307
Ω2 are filled with regular quadrilaterals while Ω3 is filled with unstructured triangles. The global hybrid mesh has
308
1888 triangles, 1204 quadrilaterals, 5044 conforming interfaces, 320 non-conforming hybrid interfaces and 120 non-
309
conforming interfaces between quadrilaterals. The average edge length of this mesh is of 0.16λ. The triangular mesh is
310
obtained by replacing the quadrilateral elements of the hybrid mesh by unstructured triangles. The resulting triangular
311
mesh is conformal and consists of 4376 triangles and 6626 interfaces with an average edge length of 0.18λ. For the
312
hybrid DGTD-P p Qq method, the interpolation order, q, is a pair q = (q1 , q2 ), where q1 and q2 are the interpolation
313
orders used in the smallest and biggest quadrilaterals. The final simulation time has been set to T = 33.33 ns which
314
corresponds to 10 periods of the incident wave. The results are listed in Table 1 in terms of accuracy and computational
315
cost to reach the final simulation time for different p and q. Inspecting these results, we observe that for a fixed p,
316
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
317
the triangular method. On Figure 9 we compare the x-wise 1D distributions of the real part of the discrete Fourier
318
transform of Ez and Hy components. We observe from this figure that, the main error of DG solutions arises at the
319
material interfaces where the solution loses smoothness and the results confirm the ability of the hybrid DG method to
320
accurately model problems with discontinuous coefficients and solutions. Finally, contour plots of the Ez component
321
at time t = T for numerical simulations performed with the DG method on the hybrid and the triangular meshes, are
322
illustrated in Figure 10 for p = 2 and q = (2, 2). We observe very good agreement between the solutions obtained
323
by both DG methods and the exact solution even with relatively low grid resolutions of six points per λ. The hybrid
324
method shows slightly better results than the triangular method.
Ac ce pt e
d
M
an
us
cr
ip t
305
1.6
1.6
0.8
0.8
0
0
-0.8
-0.8
-1.6
-1.6 -1.6
-0.8
0
0.8
1.6
-1.6
-0.8
0
0.8
1.6
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.
16
Page 16 of 28
Table 1: Comparison between the DGTD-P p Qq 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 (q1 , q2 ) where q1 and q2 are the interpolation orders used in the smallest and biggest quadrilaterals, respectively. We give the maximum ∆t allowed by the stability condition, the L2 -errors at time t = 33.33 ns, the #DOF and the
L2 -error
#DOF
CPU (s)
DGTD-P2
3.648E-02
2.394E-01
26256
201
DGTD-P3
2.132E-02
3.931E-02
43760
581
DGTD-P4
1.684E-02
6.532E-03
65640
1108
DGTD-P5
1.403E-02
1.083E-03
91896
1886
DGTD-P6
1.038E-02
1.785E-04
122528
3488
DGTD-P2 Q(2,2)
3.648E-02
1.692E-01
22164
165
DGTD-P3 Q(3,2)
2.346E-02
2.398E-02
32768
379
DGTD-P4 Q(3,2)
1.841E-02
5.105E-03
42208
592
DGTD-P5 Q(4,3)
1.523E-02
5.784E-04
62836
1162
DGTD-P6 Q(4,3)
1.216E-02
1.484E-04
76052
1671
cr
an
Method
ip t
∆t (ns)
us
computing time (CPU) in seconds to achieve the final time period.
Exact Ez DG-P3 DG-P3Q(3,2)
1
M
8
1.5
Exact Hy DG-P3 DG-P3Q(3,2)
6 4
-0.5
-1
-1.5 -1.8
d -0.6
x
0.6
field Hy
0
Ac ce pt e
field Ez
0.5
2 0 -2 0.7
zoom-in
1.2
zoom-in
-4 -6 -0.1 -0.7
-8 -1.8
1.8
0.9 0.6
-0.6
-0.6
x
0.6
0.7
1.8
Figure 9: Plots of the Ez (left) and Hy (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-P3 Q(3,2) (square) and triangular DGTD-P3 (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.
325
5.1.3. Standing wave in a wedge-shaped PEC resonator
326
As an example with geometric singularities, we consider the propagation of a standing wave of frequency f =
327
1293 MHz and wavelength λ = 0.232 m inside a domain bounded by the curves y = tan(7π/4)x, x2 + y2 = 1/4 and the
328
x-axis. The non-convex nature of the outer boundary renders some components of the solution highly singular, making
329
it a challenging problem. We compare approximate solutions obtained using the DGTD-P p method formulated on 17
Page 17 of 28
1
0.5
0.5
0.5
0
0
0
-0.5
-0.5
-0.5
-1 Y -1 Z X
-0.5 0 0.5 Ez component - exact solution
1
-1 Y -1 Z X
-0.5 0 0.5 Ez component on hybrid mesh for p=2, q=(2,2)
1
-1 Y -1 Z X
ip t
1
-0.5 0 0.5 Ez component on triangular mesh for p=2
cr
1
1
Figure 10: Contour lines of the Ez component at time t = 33.33 ns for the exact solution (left) and for solutions resulting from the hybrid DGTD-
us
P2 Q(2,2) (center) and triangular DGTD-P2 (right) methods. Zoomed-in region to show more details. In all the plots, we use 13 equally spaced contours from Ez = −1.5166 to Ez = 1.5299.
triangular meshes with those resulting from the DGTD-P p Qq method formulated on hybrid non-conforming meshes.
331
The triangular mesh consists of 6280 elements and 9531 interfaces with an average edge length of 0.057λ. The hybrid
332
mesh contains 1745 triangles, 1820 quadrilaterals, 6228 conforming interfaces and 280 non-conforming interfaces.
333
The global hybrid mesh has an average edge length of 0.082λ. Moreover, we can see that the hybrid mesh is locally
334
refined in the vicinity of the re-entrant corner located at the origin. The resulting meshes are depicted in Figure 11.
335
The simulation time has been set to T = 38.5 ns which corresponds to 50 periods. Results are summarized in Table 2
336
in terms of accuracy and computational costs. Here, the hybrid DG method produces numerical solutions that are one
337
order of magnitude more accurate than those delivered by the triangular DG method. Figure 12 shows the contour lines
338
of the H x component at t = T for numerical solutions obtained using the DGTD-P4 and DGTD-P2 Q4 method. One
339
can clearly observe a good qualitative agreement between the exact solution and the hybrid DGTD-P2 Q4 solutions.
340
The triangular DGTD-P4 solution produces some oscillations near the re-entrant corner due to bad resolution. This
341
is confirmed in Figure 13 where we plot the x-wise 1D distributions of the H x and Hy components along the line
342
y = 5x + 0.025 passing nearby the corner singularity. Inspecting Table 2 and Figure 13, we observe that, for the
343
DGTD-P p method, the error decreases slowly when increasing the order from p = 2 to p = 4. In contrast, the hybrid
344
method offers faster convergence when increasing the order in the quadrilateral elements from q = 2 to q = 4. More
345
specifically, the global error depends on the global regularity of the exact solution which motivates the use of lower
346
polynomial order together with mesh refinement near singular points in order to recover standard convergence rates.
347
5.2. Three-dimensional examples
348
5.2.1. Eigenmode in a spherical cavity
Ac ce pt e
d
M
an
330
349
As a first verification of the general 3D framework, let us consider the propagation of the (0, 1, 1)-mode inside
350
a PEC spherical cavity of unit radius. The resonant frequency is 131 MHz and the wavelength is λ = 2.29 m. The
351
exact time-domain solution of this problem can be found in [3]. We compare the accuracy and efficiency of the DG 18
Page 18 of 28
0.5
0
0
-0.5 -0.5
-0.5 -0.5
0.5
0
0.5
us
0
cr
ip t
0.5
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.
an
Table 2: Comparison between the hybrid DGTD-P p Qq method and the triangular DGTD-P p method. We give the maximum ∆t allowed by the
L2 -error
#DOF
CPU (s)
DGTD-P2
1.570E-02
5.482E-02
37680
1604
DGTD-P4
7.246E-03
2.381E-02
94200
9832
DGTD-P2 Q2
1.724E-02
9.875E-03
26850
1211
DGTD-P2 Q4
8.701E-03
1.043E-03
55970
6539
Method
-0.18 -0.18
0.18
Ac ce pt e
0.18
d
∆t (ns)
M
stability condition, the L2 -errors at time t = 38.5 ns, the #DOF and the computing time (CPU) in seconds to achieve the final time period.
0.18
Hx component - exact solution
-0.18 -0.18
0.18
0.18
Hx component on hybrid mesh for p=2, q=4
-0.18 -0.18
0.18
Hx component on triangular mesh for p=4
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 DGTDP2 Q4 (center) and triangular DGTD-P4 (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.
352
method using hexahedral, tetrahedral and hybrid meshes. First, a fully hexahedral mesh consisting of 6413 hexahedra
353
and 19602 quadrilateral faces is generated. Then, the fully tetrahedral mesh is obtained by considering the hexahedral
354
mesh and splitting each hexahedron to six tetrahedra. The resulting tetrahedral mesh has 38478 tetrahedra and 77682 19
Page 19 of 28
0.25
0.4
Exact Hx DGTD-P2Q4 DGTD-P4
0.3
Exact Hy DGTD-P2Q4 DGTD-P4
0.2
0.2
zoom-in 0.11
0.15 0.1 -0.004
field Hy
field Hx
0.1 0
0
0.1 0.05
ip t
-0.1 0
-0.2
-0.05
-0.3 -0.4
-0.1 -0.05
0
0.05
0.1
-0.1
-0.05
0
0.05
0.1
cr
-0.1
x
x
Figure 13: Plots of the H x (left) and Hy (right) components along the line y = 5x + 0.025 at time t = 38.5 ns for the wedge-shaped cavity. We
us
compare the exact solution (solid line) with solutions obtained using the DGTD-P2 Q4 (dashed line) and DGTD-P4 (dotted line) methods.
triangular faces. Finally, the hybrid mesh contains 13068 tetrahedra, 216 regular hexahedra, 26136 triangular faces,
356
540 quadrilateral faces and 2916 hybrid non-conforming interfaces. These meshes are illustrated in Figure 14. The
357
average edge lengths of the hexahedral, tetrahedral and hybrid meshes are respectively 0.0387λ, 0.0463λ and 0.0676λ
358
which correspond to 26, 22 and 15 points per wavelength. We report on results obtained by the hexahedral DGTD-Qq ,
359
tetrahedral DGTD-P p and hybrid DGTD-P p Qq methods with interpolation orders p and q ranging from p = 1 up to
360
p = 5 with p = q. The simulation time is fixed to T = 76 ns which corresponds to 10 periods. Table 3 summarizes the
361
values of time steps used in the simulation. In Figure 15 (left) we show the L2 error as a function of the interpolation
362
order. Inspecting these results we observe several things. First, the gain in accuracy is notable when increasing the
363
interpolation order for all methods. Here, the exponential decay of the error as exp(−2p) is observed for all cases.
364
Second, the tetrahedral method is 3.9 to 7.9 times more accurate than the hexahedral method when increasing the
365
interpolation order from p = 1 to p = 5. The results reveal that the accuracy of the hybrid method is almost identical
366
to that of the tetrahedral method. Figure 15 (center) shows the ratio between the #DOF required by the hexahedral,
367
tetrahedral and hybrid methods to achieve the final simulation time. The results indicate that the hybrid method is able
368
to reduce the #DOF for given p and q. For instance, with p = q = 5, the L2 error is 7.1 × 10−5 for the hybrid mesh,
369
4.7 × 10−4 for the hexahedral method and 6.0 × 10−5 for the tetrahedral mesh. With the hybrid mesh, the #DOF is
370
reduced by a factor of 1.8 compared to hexahedral mesh and by factor of 2.8 compared to tetrahedral mesh. Note that,
371
the higher #DOF required by the tetrahedral mesh stems from the fact that the number of elements in the tetrahedral
372
mesh is 6 times higher than those in the hexahedral mesh. Finally, the comparison in terms of computing times is
373
displayed in Figure 15 (right). It can be observed that the reduction in the #DOF is translated into a corresponding
374
reduction in computational cost. Here, the hybrid method allows for a reduction in computing time by a factor ranging
375
from 2.4 to 8.0 compared to hexahedral method and by a factor of 3.5 compared to tetrahedral method. The lower
376
computing times of the hybrid method partially stems from its larger time steps employed in hybrid mesh compared
377
to those used in tetrahedral and hexahedral meshes.
Ac ce pt e
d
M
an
355
20
Page 20 of 28
ip t cr
Y X
Y Z
X
Y
Z
Z
X
us
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.
an
Table 3: The time steps allowed by the stability conditions of the DG method using hexahedral, tetrahedral and hybrid meshes for different
p=q=2
Hexahedra
1.970E-02
9.850E-03
Tetrahedra
2.357E-02
1.178E-02
Hybrid
3.178E-02
1.589E-02 3
HEX TET HYB
Ratio of the DOF
10-2 10-3
10
-4
10
-5
p=q=5
5.910E-03
4.728E-03
3.546E-03
7.072E-03
5.658E-03
4.243E-03
9.535E-03
7.628E-03
5.721E-03
8
HEX/TET HEX/HYB
Ac ce pt e
L2 error
10-1
p=q=4
2
Ratio of computing times
100
p=q=3
d
p=q=1
M
interpolation orders.
TET/HYB
1
2
3
Interpolation order
4
5
1
HEX/TET HEX/HYB
6
TET/HYB
5 4 3 2 1
0
1
7
0
2
3
4
Interpolation order
5
1
2
3
4
5
Interpolation order
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.
378
5.2.2. Exposure of human head tissues to a localized source radiation
379
As a final example, we consider a more challenging problem which is concerned with the propagation of an EM
380
wave emitted by a localized source in a realistic geometrical model of head tissues. Starting from magnetic resonance
381
images of the Visible Human 2.0 project [46], head tissues are segmented and the interfaces of a selected number of
382
tissues are triangulated. Different strategies can be used in order to obtain a smooth and accurate segmentation of 21
Page 21 of 28
head tissues and interface triangulations as well. The strategy adopted in this example consists in using a variant of
384
Chew’s algorithm [9], based on Delaunay triangulation restricted to the interface, which allows to control the size and
385
aspect ratio of interface triangles. Example of triangulations of the skin, skull and brain are shown in Figure 16. Then,
386
these triangulated surfaces are used as inputs for the generation of volume meshes of the tissues. Finally, the GHS3D
387
tetrahedral mesh generator [33] is used to mesh the volume domains between the various interfaces. The exterior of
388
the head must also be meshed, up to a certain distance from the skin. The computational domain is here artificially
389
bounded by a surface, denoted by Γa , on which the absorbing boundary condition is imposed. Overall, the constructed
390
geometrical model considered here consists of four tissues (skin, skull, CSF - Cerebro Spinal Fluid and brain). The
391
characteristics of the tissues are summarized in Table 4 where the values of the relative permittivity εr , and the
392
conductivity σ, correspond to a frequency F = 1800 MHz and have been obtained from a special purpose online data
393
base. For all tissues, the relative permeability, µr , is set to 1. Given this tetrahedral mesh of the head tissues, the hybrid
394
tetrahedral-hexahedral mesh is built in the following way (see the left panel of Figure 17): (a) the artificial boundary
395
Γa is defined as a plane-parallel surface; (b) a second plane-parallel surface, denoted by Γ0 , lying close to the head
396
and completely enclosing it, is introduced; (c) the volume space between Γa and Γ0 (vacuum space) is discretized
397
using a regular hexahedral mesh; (d) the volume space between Γ0 and the skin is discretized with unstructured
398
tetrahedral elements. The global hybrid mesh consists of 197691 tetrahedra, 7000 hexahedra, 394397 triangular
399
interfaces, 21900 quadrilateral interfaces and 5124 hybrid non-conforming interfaces. The minimum, maximum and
400
average edge lengths of the hybrid mesh are equal to 1.12 mm, 43.14 mm and 17.26 mm, respectively. To compare the
401
performance of the DG methods on hybrid and tetrahedral meshes, a fully tetrahedral mesh is generated by considering
402
the hybrid mesh and replacing the hexahedral part by unstructured tetrahedral elements. We stress that the resulting
403
tetrahedral mesh is conforming i.e., the same triangulated surface Γ0 is used for the generation of tetrahedral elements
404
inside and outside Γ0 , as shown in Figure 17 (right). The global tetrahedral mesh consists of 350581 tetrahedra and
405
709395 interfaces. The minimum, maximum and average edge lengths of the tetrahedral mesh are equal to 0.49 mm,
406
44.32 mm and 13.98 mm, respectively. Finally, a dipolar type source is localized near the right ear of the head
407
yielding a current of the form Jz (x, t) = δ(x − xd ) f (t), where f (t) is a sinusoidally varying temporal signal and xd is
408
the localization point of the source. The physical simulation time has been fixed to T = 2.78 ns which corresponds to
409
five periods of the temporal signal. A discrete Fourier transform of the components of the electric field is computed
410
during the last period of the simulation. For this problem, we consider using the DGTD-P2 method on the tetrahedral
411
mesh and the DGTD-P2 Q2 and DGTD-P2 Q1 methods on the hybrid mesh. The time steps used in the simulations
412
are ∆t = 2.728E-04 ns for the DGTD-P2 method, ∆t = 6.655E-04 ns for the DGTD-P2 Q2 method and ∆t = 1.012E-
413
03 ns for the DGTD-P2 Q1 method. These values correspond to the maximum ∆t allowed by the stability conditions
414
of the DG methods. Contour lines of the modulus of the electric field, |E|, on selected planes of the skin for the
415
approximate solutions resulting from the hybrid DGTD-P2 Q2 and tetrahedral DGTD-P2 methods are visualized in
416
Figure 18. Visually both methods provide almost the same solutions except slight amplitude differences around the
417
right ear of the head. For the sake of completeness, we compare in Figure 19 the time evolution of the Ez and Hy
Ac ce pt e
d
M
an
us
cr
ip t
383
22
Page 22 of 28
components at a selected point in the center of the head. The comparisons show a visually perfect match between all
419
solutions and no numerical artifacts due to the non-conformity of the mesh is observed. In Figure 19 we also plot
420
the differences between the DGTD-P2 solution and the hybrid solutions, showing that the errors are relatively very
421
small (less than 0.1% for DGTD-P2 Q1 and less than 0.01% for DGTD-P2 Q2 ). Finally, we give some informations on
422
the simulation times required by both schemes. On a workstation equipped with an Intel Xeon CPU E5-2650 @2.00
423
GHz and 64 GB of RAM, the DGTD-P2 method requires 27 h 47 min for a total of 10183 time steps, the DGTD-P2 Q2
424
method requires 7 h 35 min for a total of 4177 time steps, and the DGTD-P2 Q1 method requires 4 h 17 min for a total
425
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
426
DG method.
M
an
us
cr
ip t
418
Z
Z
YX
Z Y
X
d
YX
Ac ce pt e
Figure 16: Surface mesh of the skin (left), the skull (center) and the brain (right).
Table 4: Electromagnetic characteristics of the selected head tissues at frequency f = 1800 MHz.
Tissue Skin Skull CSF Brain
427
εr
σ (S/m)
wavelength (mm)
38.87
1.18
26.73
15.56
0.43
42.25
67.20
2.92
20.33
43.55
1.15
25.26
6. Concluding remarks
428
We have reported on our efforts to design a high-order accurate DG method for solving the time-domain Maxwell
429
equations on non-conforming hybrid meshes. The proposed method combines a central numerical flux with a fourth-
430
order leap-frog time integration. Within each mesh element, the EM field components are approximated using ar-
431
bitrary high-order nodal polynomials. The geometrical flexibility is achieved by using unstructured triangular or 23
Page 23 of 28
ip t cr us
an
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.
tetrahedral meshes around complex geometric details and computational efficiency is then improved by using quadri-
433
laterals or hexahedra to fill the remainder of the domain. The resulting DGTD method conserves a discrete form of
434
the EM energy, and is stable under some CFL-like conditions. The key issues for the practical implementation of
435
the proposed method are the generation of the hybrid meshes and the computation of flux matrices at the interfaces
436
between hexahedral and tetrahedral faces. First, the mesh strategy adopted in this paper relies on the decomposition
437
of the computational domain into non-overlapping subdomains, each one is meshed independently using the appro-
438
priate element type. Then, the hybrid interfaces are found using the method of separating axis and Boolean operations
439
over polygons. This strategy demonstrates the benefits of relaxing the conformity requirement in the discretization
440
process in view of facilitating the construction of meshes for complex propagation configurations, and improving the
441
overall efficiency of the simulations. It is noteworthy that, most, if not all, of commercial or non-commercial mesh-
442
generating softwares currently produce only conforming hybrid meshes by adding transition layers of pyramids or
443
prisms. Second, the flux matrices are calculated using cubature rules over polygons. For a given accuracy, these rules
444
are able to save more than 65% of the cubature points compared to other approaches. The proposed DGTD method is
445
implemented in two and three space dimensions and tests on different wave propagation problems involving complex
446
realistic geometries and heterogeneous media are presented. The spectral convergence for p-refinement method and
447
the optimal convergence rate for h-refinement method have been verified. Results also show excellent agreement with
448
exact solutions and no numerical artifacts caused by the non-conformity of the mesh are observed. Compared to DG
449
methods on single element-type meshes, the DG method on hybrid meshes leads to clear reduction in computing time
450
which is highly dependent on the interpolation order. In 2D, it is found that the hybrid mesh is able to produce at least
451
a similar accuracy and allows the DG method to save up to 55% of the computing time and the number of degrees of
Ac ce pt e
d
M
432
24
Page 24 of 28
0.514
0.481
0.481
0.447
0.447
0.413
0.413
0.38
0.38
0.346
0.346
0.313
0.313
0.279
0.279
0.245
0.245
0.212
0.212
ip t
0.514
0.178
0.178
0.144
0.144
0.111
0.111
0.0772
0.0772
0.0436
0.504 0.471 0.438 0.406 0.373 0.307 0.274 0.241 0.209
0.0772 0.0443 0.0115
cr
Y Z X
|E| - DGTD-P2Q2
0.504 0.471 0.438 0.406 0.373 0.34 0.307 0.274 0.241 0.209 0.176 0.143 0.11 0.0772 0.0443 0.0115 |E| - DGTD-P2
d
Y Z X
M
0.176 0.11
|E| - DGTD-P2
an
0.34
0.143
0.01
Z Y X
|E| - DGTD-P2Q2
us
Z Y X
0.0436
0.01
Ac ce pt e
Figure 18: Contour lines of the modulus of the electric field, |E|, in selected cut planes of the head, using the hybrid DGTD-P2 Q2 method (left) and the tetrahedral DGTD-P2 method (right).
452
freedom. In 3D, the gain becomes more significant, as the computational cost is reduced by up to 85%. A comparison
453
between a fully hexahedral mesh and a fully tetrahedral mesh (obtained by converting each hexahedron to six tetrahe-
454
dra) reveals that the tetrahedral method is always more accurate and more efficient than hexahedral method, excepting
455
for linear interpolation orders. The parallelization of the proposed hybrid DG method and the comparison with DG
456
methods on conforming hybrid meshes are left for future study.
457
Acknowledgements
458
The author would like to thank Dr. Seyed Mousavi of the University of California at Davis, USA, for providing
459
the quadrature points over polygons to calculate flux matrices at hybrid interfaces. We also thank Dr. Nicolas Roquet
460
for several discussions concerning the intersection tests for hybrid mesh generation. Finally, the author acknowledge
461
partial support from the French National Agency for Research (ANR-10-JCJC-0905).
25
Page 25 of 28
0.05
0.1 zoom-in
0.043
0.04
0.08
0.03
0.06 0.04
0.01 2.57
2.65
0 -0.01
0.02 0 -0.02
-0.02
DGTD-P2
-0.03
DGTD-P2Q2
-0.04
DGTD-P2
-0.06
DGTD-P2Q2
DGTD-P2Q1
DGTD-P2Q1 -0.08
cr
-0.04
-4
-3
10
10
(P2 - P2Q1)
0
0
(P2 - P2Q2) 0
0.5
(P2 - P2Q2)
1 1.5 2 time (nanoseconds)
-10
2.5
-4
0
0.5
an
-3
us
error
error
(P2 - P2Q1)
-10
ip t
0.037
field Ez
field Hy
0.02
1 1.5 2 time (nanoseconds)
2.5
Figure 19: The top panels show the time evolution of Hy (left) and Ez (right) components at a selected point in the head. These solutions are obtained
M
using the fully tetrahedral DGTD-P2 method (solid line) and the hybrid DGTD-P2 Q2 (dashed line) and DGTD-P2 Q1 (dotted line) methods. In the bottom panels we plot the differences between the tetrahedral and hybrid DG solutions denoted as P2 − P2 Q2 and P2 − P2 Q1 .
464 465 466
d
463
References
[1] E. Abenius, U. Andersson, F. Edelvik, L. Eriksson, G. Ledfelt, Hybrid time domain solvers for the maxwell equations in 2d, Int. J. Numer. Meth. Eng. 53 (9) (2002) 2185–2199.
Ac ce pt e
462
[2] T. Apel, N. Düvelmeyer, Transformation of hexaedral finite element meshes into tetrahedral meshes according to quality criteria, Computing 71 (4) (2003) 293–304.
467
[3] C. A. Balanis, Advanced Engineering Electromagnetics, John Wiley & Sons, New York, 1989.
468
[4] K. Barrett, Multilinear jacobians for isoparametric planar elements, Finite Elem. Anal. Des. 40 (8) (2004) 821–853.
469
[5] M. Bergot, M. Duruflé, Higher-order discontinuous Galerkin method for pyramidal elements using orthogonal bases, Numer. Meth. Part. Diff.
470 471 472
Eq. 29 (1) (2013) 144–169.
[6] L. Branets, G. Carey, Extension of a mesh quality metric for elements with a curved boundary edge or surface, J. Comput. Inf. Sci. Eng. 5 (4) (2005) 302–308.
473
[7] K. Busch, M. König, J. Niegemann, Discontinuous Galerkin methods in nanophotonics, Laser & Photon. Rev. 5 (6) (2011) 773–809.
474
[8] G. Carey, Hexing the Tet, Communications in Numerical Methods in Engineering 18 (3) (2002) 223–227.
475
[9] L. Chew, Guaranteed-quality mesh generation for curved surfaces, in: 9th Annual ACM Symposium Computational Geometry, ACM Press,
476 477 478 479 480 481
1993, pp. 274–280.
[10] B. Cockburn, G. Karniadakis, C. Shu (eds.), Discontinuous Galerkin Methods. Theory, Computation and Applications, vol. 11 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin, 2000. [11] G. Cohen, X. Ferrieres, S. Pernet, A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell’s equations in time domain, J. Comput. Phys. 217 (2) (2006) 340–363. [12] R. Cools, Monomial cubature rules since “Stroud”: a compilation–part 2, J. Comput. Appl. Math. 112 (1-2) (1999) 21–27.
482
[13] R. Cools, An encyclopaedia of cubature formulas, J. Complex. 19 (3) (2003) 445–453, online database (http://nines.cs.kuleuven.be/ecf/).
483
[14] R. Cools, P. Rabinowitz, Monomial cubature rules since “Stroud”: a compilation, J. Comput. Appl. Math. 48 (3) (1993) 309–326.
26
Page 26 of 28
484 485 486 487
[15] S. Dosopoulos, B. Zhao, J.-F. Lee, Non-conformal and parallel discontinuous galerkin time domain method for maxwell’s equations: EM analysis of IC packages, J. Comput. Phys. 238 (0) (2013) 48–70. [16] D. Dunavant, Economical symmetrical quadrature rules for complete polynomials over a square domain, Int. J. Numer. Meth. Eng. 21 (10) (1985) 1777–1784.
488
[17] D. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle, Int. J. Numer. Meth. Eng. 21 (6) (1985) 1129–1148.
489
[18] D. Dunavant, Efficient symmetrical cubature rules for complete polynomials of high degree over the unit cube, Int. J. Numer. Meth. Eng.
491 492
23 (3) (1986) 397–407.
ip t
490
[19] C. Durochat, S. Lanteri, DGTD method on hybrid meshes for time domain electromagnetics, in: URSI International Symposium on Electromagnetic Theory (EMTS), 2010, pp. 992–995.
[20] C. Durochat, S. Lanteri, C. Scheid, High order non-conforming multi-element discontinuous Galerkin method for time-domain electromagnetics, in: International Conference on Electromagnetics in Advanced Applications (ICEAA), Cape Town, South Africa, 2012, pp. 379–382.
cr
493 494
[21] D. Eberly, 3D Game Engine Design, 2nd ed., Morgan Kaufmann Publishers, 2007.
496
[22] D. Eberly, Game Physics, 2nd ed., Morgan Kaufmann Publishers, 2010.
497
[23] F. Edelvik, G. Ledfelt, Explicit hybrid time domain solver for the Maxwell equations in 3D, J. Sci. Comput. 15 (1) (2000) 61–78.
498
[24] H. Fahs, Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic
503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521
an
502
[25] H. Fahs, High-order leap-frog based discontinuous Galerkin method for the time-domain Maxwell equations on non-conforming simplicial meshes, Numer. Math. Theor. Meth. Appl. 2 (3) (2009) 275–300.
[26] H. Fahs, Improving accuracy of high-order discontinuous Galerkin method for time-domain electromagnetics on curvilinear domains, Int. J. Comput. Math. 88 (10) (2011) 2124–2153.
M
501
[27] H. Fahs, Investigation on polynomial integrators for time-domain electromagnetics using a high-order discontinuous Galerkin method, Appl. Math. Model. 36 (11) (2012) 5466–5481.
[28] H. Fahs, L. Fezoui, S. Lanteri, F. Rapetti, Preliminary investigation of a non-conforming discontinuous Galerkin method for solving the
d
500
wave propagation, Int. J. Numer. Anal. Model. 6 (2) (2009) 193–216.
time-domain Maxwell equations, IEEE Trans. Magn. 44 (6) (2008) 1254–1257. [29] H. Fahs, A. Hadjem, S. Lanteri, J. Wiart, M.-F. Wong, Calculation of the SAR induced in head tissues using a high order DGTD method and
Ac ce pt e
499
us
495
triangulated geometrical models, IEEE Trans. Antennas Propag. 59 (12) (2011) 4669–4678. [30] H. Fahs, S. Lanteri, A high-order non-conforming discontinuous Galerkin method for time-domain electromagnetics, J. Comput. Appl. Math. 234 (4) (2010) 1088–1096.
[31] X. Ferrieres, J.-P. Parmantier, S. Bertuol, A. R. Ruddle, Application of a hybrid finite difference/finite volume method to solve an automotive EMC problem, IEEE Trans. Electromagn. Compat. 46 (4) (2004) 624–634. [32] G. Gassner, F. Lörcher, C.-D. Munz, J. Hesthaven, Polymorphic nodal elements and their application in discontinuous Galerkin methods, J. Comput. Phys. 228 (5) (2009) 1573–1590.
[33] P.-L. George, F. Hecht, E. Saltel, Automatic mesh generator with specified boundary, Comput. Methods Appl. Mech. Engrg. 92 (3) (1991) 269–288.
[34] S. Gottschalk, M. Lin, D. Manocha, OBB-tree: a hierarchical structure for rapid interference detection, in: Proceedings of ACM Siggraph’96, ACM, New York, USA, 1996, pp. 171–180.
[35] J. Hesthaven, From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex, SIAM J. Numer. Anal. 35 (2) (1998) 655–676.
522
[36] J. Hesthaven, C. Teng, Stable spectral methods on tetrahedral elements, SIAM J. Sci. Comput. 21 (6) (1999) 2352–2380.
523
[37] J. Hesthaven, T. Warburton, Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell’s equations, J. Comput.
524 525 526
Phys. 181 (1) (2002) 186–221. [38] J. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods: Algorithms, Analysis and Applications, Springer Texts Appl. Math., Springer-Verlag, Berlin, 2008.
27
Page 27 of 28
527
[39] P. Hubbard, Constructive solid geometry for triangulated polyhedra, Tech. Rep. CS-90-07, Brown University (1990).
528
[40] P. Knabner, S. Korotov, G. Summ, Conditions for the invertibility of the isoparametric mapping for hexahedral finite elements, Finite Elem.
529
Anal. Des. 40 (2) (2003) 159–172.
530
[41] P. M. Knupp, On the invertibility of the isoparametric map, Comput. Methods Appl. Mech. Engrg. 78 (3) (1990) 313–329.
531
[42] T. C. L. Zhang, H. Liu, A set of symmetric quadrature rules on triangles and tetrahedra, J. Comp. Math. 27 (1) (2009) 89–96.
532
[43] S. Mousavi, H. H. Xiao, N. Sukumar, Generalized gaussian quadrature rules on arbitrary polygons, Int. J. Numer. Meth. Eng. 82 (1) (2010)
536 537 538 539
ip t
535
[44] S. Pernet, X. Ferrieres, G. Cohen, High spatial order finite element method to solve Maxwell’s equations in time domain, IEEE Trans. Antennas Propag. 53 (9) (2005) 2889–2899.
[45] S. Piperno, M. Remaki, L. Fezoui, A nondiffusive finite volume scheme for the three-dimensional Maxwell’s equations on unstructured
cr
534
99–113.
meshes, SIAM J. Numer. Anal. 39 (6) (2002) 2089–2108.
[46] P. Ratiu, B. Hillen, J. Glaser, D. Jenkins, Medicine Meets Virtual Reality 11 - NextMed: Health Horizon, vol. 11, chap. Visible Human 2.0 the next generation, IOS Press, 2003, pp. 275–281.
us
533
540
[47] M. Rivero, F. Feito, Boolean operations on general planar polygons, Computers & Graphics 24 (6) (2000) 881–896.
541
[48] T. Rylander, A. Bondeson, Stability of explicit-implicit hybrid time-stepping schemes for Maxwell’s equations, J. Comput. Phys. 179 (2) (2002) 426–438.
an
542 543
[49] P. Schneider, D. Eberly, Geometric Tools for Computer Graphics, 1st ed., Morgan Kaufmann Publishers, 2003.
544
[50] S. Schnepp, T. Weiland, Efficient large scale electromagnetic simulations using dynamically adapted meshes with the discontinuous Galerkin
547 548
method, J. Comput. Appl. Math. 236 (18) (2012) 4909 – 4924.
[51] R. Sevilla, O. Hassan, K. Morgan, The use of hybrid meshes to improve the efficiency of a discontinuous galerkin method for the solution of maxwell’s equations, Comput. Struct.In Press.
M
545 546
[52] H. Songoro, M. Vogel, Z. Cendes, Keeping time with Maxwell’s equations, IEEE Microw. Mag. 11 (2) (2010) 42–49. [53] N. Sukumar, A. Tabarraei, Conforming polygonal finite elements, Int. J. Numer. Meth. Eng. 61 (12) (2004) 2045–2066.
550
[54] A. Taflove, S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed., Artech House, Norwood, MA, 2005.
Ac ce pt e
551
d
549
552
[55] S. Wandzura, H. Xiao, Symmetric quadrature rules on a triangle, Comput. Math. Appl. 45 (12) (2003) 1829–1840.
553
[56] R. B. Wu, T. Itoh, Hybrid finite-difference time-domain modeling of curved surfaces using tetrahedral edge elements, IEEE Trans. Antennas
554
Propag. 45 (8) (1997) 1302–1309.
28
Page 28 of 28