# Electromagnetic Theory with Discrete Exterior Calculus.

1. INTRODUCTION

Finite difference (FD) and finite element methods (FEM) have been widely applied in electromagnetics [1, 2]. And they are both based on the vectorial version of Maxwell's equations. As is well known, differential forms can be used to recast Maxwell's equations in a more succinct fashion, which completely separates metric-free and metric-dependent components [4-9]. Discrete exterior calculus (DEC), a discrete version of differential geometry, provides a numerical treatment of differential forms directly [10-12]. Compared to the popular finite difference method, DEC is based on unstructured simplicial mesh as in FEM, e.g., triangular mesh for 2-D space and tetrahedral mesh for 3-D space. Therefore, this method is adaptable to arbitrary complex structures. Moreover, as will be shown, this method also exactly preserves important structural features of Maxwell's equations as in Yee gird applied to FD, e.g., Gauss's law [nabla] x D = [rho] and E = -[nabla][phi]. Meanwhile, unlike finite integration technique (FIT) [15,16], the dual mesh is constructed by connecting circumcenters instead of barycenters. In this way, the Hodge star operators can be built as diagonal matrices. Effective dual volume can be applied to incorporate the material information to include inhomogeneity as shown in .

Shown in , we have been able to solve eigenvalue problems, excitation problems, and scattering problems with DEC method. It should be noted that the only limitation of our approach is that without Delaunay triangulation, the positivity of Laplacian cannot be guaranteed [13,14], which will lead to instability of time domain analysis. However it is still of interest to investigate if a self-consistent and self-contained discrete electromagnetic theory can be developed on simplicial mesh using DEC, which has been accomplished based on Yee grid [3,4]. In this paper, we show that many electromagnetic theorems still exist in this discrete electromagnetic theory with DEC, such as Huygens' principle, reciprocity theorem, Poynting's theorem and uniqueness theorem. Therefore, all electromagnetic phenomena can be predicted and explained with this discrete electromagnetic theory with a appropriate discretization length scale. Besides, due to the charge conservation property, as will be shown, the numerical simulation of Maxwell's equations with DEC will not lead to spurious solutions [3,17].

To develop this discrete electromagnetic theory with DEC, an appropriate discrete treatment of wedge product [11,18] is crucial. The wedge product operator [conjunction] is used to construct higher degree forms from lower degree ones. Correspondingly, a discrete wedge product operator should be able to obtain higher degree cochains from lower degree ones. For example, with a known primal 1-cochain E and dual 2-cochain D, an integral of the wedge product of these two cochains should represent stored electric energy density E * D. It should be noted that a wedge product between two primal or two dual cochains is discussed in . In this paper, a proper treatment of wedge product between a primal fc-cochain and a dual (n - k)-cochain is given. With this discrete wedge product, the discrete electric energy density E x D, magnetic energy density H x B, and Poynting's vector E x H can all be constructed. Moreover, the convolution of Green's function with current or charge fields can also be understood as a wedge product. Therefore, Huygens' principle for both scalar and electromagnetic waves can be derived in this discrete world. And it is straightforward to show that many electromagnetic theorems, like Poynting's theorem and reciprocity theorem are also preserved.

In Section 2, discrete calculus based on DEC is introduced, including discrete gradient, curl, and divergence for both primal and dual cochains. In Section 3, we define the discrete wedge product between primal k-cochain and dual (n - k)-cochain in n-D space. In the next section, the discrete Green's functions and Huygens' principle for scalar waves and vector waves are formulated and derived. In addition, discrete Poynting's theorem and reciprocity theorem are also derived.

2. DISCRETE CALCULUS WITH DEC

In DEC, continuous scalar or vector field are represented with cochains, which are column vectors. Then the derivative operators acting on them are replaced by co-boundary or discrete exterior derivative operators, which are highly sparse matrices. As mentioned, DEC is based on arbitrary simplicial mesh, which is called the primal mesh, and the dual mesh is constructed by connecting nearest circumcenters. Cochains defined on the primal and dual mesh is called primal cochains and dual cochains, respectively. Due to the existence of boundaries, the dual mesh is truncated. Therefore, discrete calculus for dual cochains will give rise to boundary terms.

2.1. Discrete Gradient, Curl and Divergence for Primal Cochains

With a scalar function [phi](r), a gradient can be defined. With DEC, a scalar function corresponds to a 0-cochain [phi], which lives on [N.sub.0] vertices {[mathematical expression not reproducible]}, and it can be written as a [N.sub.0] x 1 column vector:

[mathematical expression not reproducible]. (1)

If we view gradient of function [phi](r), [nabla][phi](r), as a 1-cochain e, then the value of this 1-cochain along the i-th primal edge [l.sub.i] is defined as:

[mathematical expression not reproducible] (2)

where [d.sup.(0).sub.i,j] is the (i,j) element of discrete exterior derivative operator [[bar.d].sup.(0)]. More specifically,

[mathematical expression not reproducible], (3)

where the superscript "(0)" means that it only operates on 0-cochains. In fact, the rightmost of Equation (2) is just the i-th element of 1-cochain [[bar.d].sup.(0)][phi], which means that the gradient relation e = [[bar.d].sup.(0)][phi] is naturally and exactly preserved.

2.1.2. Discrete Curl and Stokes' Theorem

A curl of a vector field, i.e., E(r), can measure its circulation or rotation about a point. In both 2-D and 3-D spaces, we view this vector field as a 1-cochain E, because only a path integral of it gives rise to a sensible physical meaning. Therefore, its curl [nabla] x E(r) corresponds to a 2-cochain b in both 2-D and 3-D spaces.

In a 2-D or 3-D simplicial mesh, if there are [N.sub.1] edges in total, 1- cochain E can be defined as:

[mathematical expression not reproducible]. (4)

Here, [[??].sub.i] is the unit vector along the i-th primal edge [l.sub.i].

As shown in Figure 1, the value of 2-cochain b on i-th primal face [S.sub.i] can be derived from Stokes' theorem as:

[mathematical expression not reproducible]. (5)

where [partial derivative][S.sub.i] refers to the three primal edges around [S.sub.i]; [[??].sub.i,j] is a unit vector surrounding [S.sub.i] clockwise and along [l.sub.j]; [d.sup.(1).sub.i,j] is the (i,j) element of discrete exterior derivative operator [[bar.d].sup.(1).sub.i,j].

[mathematical expression not reproducible], (6)

It can be seen that the rightmost of Equation (5) is the i-th element of [[bar.d].sup.(1)]E. Then the relation b = [[bar.d].sup.(1)]E holds exactly in both 2-D and 3-D.

2.1.3. Discrete Divergence and Gauss' Theorem

With a vector field D(r), its divergence can be defined in both 2-D and 3-D space. In 2-D, D(r) is viewed as a primal 1-cochain ([dagger]) and its divergence [[nabla].sub.s] x D(r) corresponds to a primal 2-cochain [sigma]. While in 3-D, D(r) is viewed as a primal 2-cochain with divergence [nabla] x D(r) corresponding to a primal 3-cochain [rho].

For 2-D triangular mesh, if there are [N.sub.1] edges, the primal 1-cochain D is defined as:

[mathematical expression not reproducible]. (7)

Here, [mathematical expression not reproducible] is a unit vector perpendicular to the i-th primal edge [l.sub.i].

Therefore, as shown in Figure 2(a), the value of 2-cochain [sigma] on i-th primal triangle [S.sub.i] can be derived by 2-D Gauss' theorem:

[mathematical expression not reproducible], (8)

where [mathematical expression not reproducible] is a unit vector pointing outside [S.sub.i] and perpendicular to [[??].sub.j], and [omega][S.sub.i], [d.sup.(1).sub.i,j] are introduced earlier.

Therefore, the rightmost of Equation (8) is the i-th element of - [[bar.d].sup.(1)]D. Then relation [sigma] = - [[bar.d].sup.(1)]D is exact. It should be pointed out that this negative sign is introduced by the definition of cochain D, shown in Equation (7).

Similarly, in a 3-D tetrahedral mesh, if there are [N.sub.2] triangular faces. The 2-cochain D is defined as:

[mathematical expression not reproducible] (9)

where [[??].sub.i] is the unit vector normal to the i-th primal face [S.sub.i].

Therefore, as shown in Figure 2(b), the value of 3-cochain [rho] on i-th tet [V.sub.i] can be derived by 3-D Gauss' theorem:

[mathematical expression not reproducible]. (10)

Here, [partial derivative][V.sub.i] is the closed surface of [V.sub.i], which includes four triangles; [[??].sub.i,j] is a unit vector pointing outside [V.sub.i] and normal to primal face [S.sub.j]; [d.sup.(2).sub.i,j] is the (i, j) element of discrete exterior derivative operator [[bar.d].sup.(2).sub.i,j], which is defined as

[mathematical expression not reproducible]. (11)

Notice that the rightmost of Equation (10) is the i-th element of [[bar.d].sup.(2)]D. Then relation [rho] = [d.sup.(2).sub.i,j] is also exact.

2.2. Discrete Gradient, Curl and Divergence for Dual Cochains

In the previous section, we demonstrate the discrete gradient, divergence, and curl for cochains on the primal mesh. Indeed, for dual cochains, similar relations can be derived. Here, dual cochains are defined on dual cells, ([double dagger]) which are constructed by connecting circumcenters, as shown in Figure 3. And notice that, Gauss' theorem and Stokes' theorem can also be applied for complex structures, such as polygons or polyhedrons. Therefore, in the same way, exact relations can also be constructed for gradient, divergence and curl for dual cochains.

However, since the dual mesh is deduced from the primal mesh by connecting nearest circumcenters, the mesh is incomplete at the boundary. For example, in Figure 3, the dual edges [bar.[d.sub.1][e.sub.1]] and [bar.[e.sub.2][d.sub.2]] are in fact "truncated" by the primal mesh boundary. And these two dual edges are not line segment between two dual vertices, because there is no well defined dual vertices outside the primal mesh boundary. Instead, dual edges on the boundary are constructed by connecting circumcenters ([d.sub.1] and [d.sub.2]) and centers of corresponding primal edges ([e.sub.1] for [bar.[p.sub.1][p.sub.2]], and [e.sub.2] for [bar.[p.sub.2][p.sub.3]]). Moreover, the dual faces and volumes surrounded by these boundary dual edges and faces are also incomplete. For example, the dual face [mathematical expression not reproducible] ([section]) in Figure 3 is only enclosed by three dual edges (solid lines): [bar.[d.sub.1][e.sub.1]], [bar.[e.sub.2][d.sub.2]], and [bar.[d.sub.2][d.sub.1]]. The other two edges (dashed lines) [bar.[e.sub.1][p.sub.2] and [bar.[p.sub.2][e.sub.2]] are half of primal edges. Due to this incompleteness of dual mesh on the boundary, a certain kind of boundary condition needs to be imposed for discrete gradient, curl, divergence of dual cochains.

2.2.1. Discrete Gradient for Dual Cochains

In a 2-D mesh, a scalar function [psi](r) can also be treated as a dual 0- cochain [psi] with length [N.sub.2] on the dual mesh:

[mathematical expression not reproducible] (12)

where [d.sub.i] is the circumcenter of i-th primal face, also the i-th dual vertex. In 3-D meshes, the length of this 0-cochain becomes [N.sub.3], the number of tetrahedrons.

Then the gradient of [psi](r), [nabla][psi](r), corresponds to a dual 1-cochain h. The value of this dual 1-cochain along dual edges strictly inside, which are not at the boundary, can be derived in the same manner as we did for primal cochains in the previous section. In a 2-D mesh, along i-th dual edge [L.sub.i]:

[mathematical expression not reproducible] (13)

Here, the unit vector along dual edge [[laplace].sub.i] is obtained by rotating unit vector along corresponding primal edge [l.sub.i] by 90 degrees anticlockwise,

[mathematical expression not reproducible]

as illustrated in Figure 3 by the direction of dual edge [bar.[d.sub.2][d.sub.1]] and primal edge [bar.[p.sub.2][p.sub.0]]. The last equality in Equation (13) is because that in a 2-D mesh below relations hold true .

[[bar.d].sup.(0).sub.dual] = [([[bar.d].sup.(1)]).sup.T] and [[bar.d].sup.(1).sub.dual] = -[([[bar.d].sup.(0)]).sup.T]. (14)

Therefore, for cochain values defined on edges strictly inside, h = [([[bar.d].sup.(1)]).sup.T] [psi] is exactly satisfied. While in 3-D, since Equation (14) becomes 

[bar.d(0).sub.dual] = [([[bar.d].sup.(2)]).sup.T], [[bar.d].sup.(1).sub.dual] = [([[bar.d].sup.(1)]).sup.T], [[bar.d].sup.(2).sub.dual] = -[([[bar.d].sup.(0)]).sup.T], (15)

Then the equation becomes h = [([[bar.d].sup.(2)]).sup.T][psi].

However, along dual edges at the boundary, such as [bar.[d.sub.1][e.sub.1]] in Figure 3, Dirichlet or Neumann boundary conditions need to be imposed. For these two kinds of boundary conditions, the value of dual 1-cochain h on the boundary dual edges can be defined as:

[mathematical expression not reproducible], (16)

[mathematical expression not reproducible]. (17)

Here, [absolute value of (...)] means the length of a line segment.

2.2.2. Discrete Curl for Dual Cochains

For a vectorial case, if we have a vector field H(r), then it can be treated as a dual 1-cochain H on the dual mesh as:

[mathematical expression not reproducible] (18)

Then the curl of this vector field [nabla] x H(r) is represented by a dual 2- cochain d, and its value on dual faces strictly inside can be derived the same as before, which means that

d = [[bar.d].sup.(1).sub.dual] H = -[([[bar.d].sup.(0)]).sup.T]H

for a 2-D mesh, and

d = [[bar.d].sup.(1).sub.dual] H = [([[bar.d].sup.(1)]).sup.T] H

for a 3-D mesh.

However, as mentioned, these relations are incomplete at the boundary. For example, in a 2-D mesh as shown in Figure 3, the vale of d on [mathematical expression not reproducible] can be derived by Stokes' theorem as:

[mathematical expression not reproducible]. (19)

Here, the first three terms are the value of dual 2-cochain [[bar.d].sup.(1).sub.dual]H or - [([[bar.d].sup.(0)]).sup.t] H, associated with dual face [mathematical expression not reproducible], while the last two terms [mathematical expression not reproducible] need to be determined by a certain boundary conditions. For example, with perfect magnetic conductor (PMC) boundary condition, they are both zero. For a general boundary condition, they have nonzero values which need to be left. Therefore, we should have a corrected relation as below:

d = -[([[bar.d].sup.(0)]).sup.T] H + [J.sub.s]. (20)

Here, [J.sub.s] is a dual 2-cochain in this 2-D mesh, and represents the tangential component of H on the boundary. In fact, if we restrict just on the boundary, which is a 1-D discrete manifold, this [J.sub.s] cochain is viewed as a dual 1-cochain. This cochain is only defined on the boundary dual faces. For example, its value on [mathematical expression not reproducible] is

[mathematical expression not reproducible]. (21)

Here, [H.sub.t] is the average tangential value of H(r) field, and [mathematical expression not reproducible]. In practice, with different boundary conditions, [J.sub.s] may have different values . When there is no boundary, which means that the space we consider is closed, although not possible in computational electromagnetics, [J.sub.s] also disappears.

It should be noted that although a 2-D example is used to illustrate the curl operation on vector field H(r), in 3-D space, the derivation would be exactly the same except that the dual face at the boundary is no longer just associated with a vertex, e.g., [mathematical expression not reproducible]. Instead, the dual face at the boundary is associated with a primal edge, and below relation holds

d = [([[bar.d].sup.(1)]).sup.T] H + [J.sub.s]. (22)

Here, [J.sub.s] is also a dual 2-cochain, and can be viewed as a dual 1-cochain on the boundary, which is a 2-D discrete manifold now.

2.2.3. Discrete Divergence for Dual Cochains

The discrete divergence for dual cochains can also be induced similarly. Certain boundary conditions are also needed to complete this relation. To be clear, we first consider the divergence of a vector field D(r) in 2-D space. This field can be viewed as a dual 1-cochain D as:

[mathematical expression not reproducible]. (23)

Here, [[??].sub.i], the unit vector along primal edge [l.sub.i], is normal to dual edge [[laplace].sub.i], shown in Figure 4. Then, by Gauss' theorem, the discrete divergence of this vector field [[nabla].sub.s] x D can be represented as a dual 2-cochain [sigma]. Strictly inside the boundary, following relation holds

[sigma] = [[bar.d].sup.(1).sub.dual]D = -[([[bar.d].sup.(0)]).sup.T] D.

At the boundary, the value of a on dual face [mathematical expression not reproducible] in Figure 4 is written as:

[mathematical expression not reproducible] (24)

where [??] is a unit vector normal to the boundary and pointing outside [mathematical expression not reproducible]. Similarly, the first three terms in the second line of Equation (24) are the value of dual 2-cochain [[bar.d].sup.(1).sub.dual]D or - [([[bar.d].sup.(0)]).sup.T] D on this dual face [mathematical expression not reproducible], but the last two terms [mathematical expression not reproducible] are not included. Therefore, we should have relation:

[sigma] = -[([[bar.d].sup.(0)]).sup.T] D + [[lambda].sub.s]. (25)

Here, [[lambda].sub.s] is a dual 2-cochain representing the normal component of D on the boundary. If only considering the boundary, which is a discrete 1-D manifold, [[lambda].sub.s] is a dual 1-cochain. And it only contains elements defined on the boundary dual faces. For example, on dual face [mathematical expression not reproducible], the value of [[lambda].sub.s] is

[mathematical expression not reproducible] (26)

where [D.sub.n] = D(r) x [??] is the average normal component of D field. Therefore, with the value of [D.sub.n] given as a boundary condition, we can determine [[lambda].sub.s].

In a 3-D mesh, the relation is quite similar. But now the vector field D(r) should be viewed as a dual 2-cochain, and its divergence is a dual 3-cochain [rho], which represents the volume charge density including a strictly inside component and a boundary (surface) component:

[rho] = -[([[bar.d].sup.(0)]).sub.T] D + [[sigma].sub.s]. (27)

Here, [[sigma].sub.s] is a dual 3-cochain, and it also only contains elements on the boundary. As shown in Figure 5, the value of [[sigma].sub.s] on [mathematical expression not reproducible] is

[mathematical expression not reproducible] (28)

where [mathematical expression not reproducible] refers to the area of [mathematical expression not reproducible]. Since the boundary surface [partial derivative][OMEGA] can be viewed as a discrete 2-D manifold, cochain as is a dual 2-cochain on this manifold. If the normal component of D(r) field [D.sub.n] at the surface is given as a boundary condition, as is determined.

2.3. Discrete Laplacian Operator

With discrete divergence and gradient defined for both primal and dual cochains, the Laplacian operator [nabla] * [nabla] can then be well defined for a scalar field [phi](r). For simplicity, we treat this scalar field as a primal 0-cochain [phi]. Then the first gradient operator take the primal 0-cochain to primal 1-cochain [[bar.d].sup.(0)] [phi]. Then, in 3-D space, for next divergence operation to function properly, a Hodge star operator is needed to map this primal 1-cochain to a dual 2-cochain:

[[bar.d].sup.(0)][phi] [right arrow] [*.sup.(1)][bar.d][phi]

where the superscript of Hodge star operator means the degree of cochains it acts on. Then as discussed earlier in Equation (27), the divergence of this dual 2-cochain will introduce a boundary term:

[nabla] x [nabla][phi](r) [??] -[[([[bar.d].sup.(0)]).sup.T] [*.sup.(1)][[bar.d].sup.(0)]] [phi] + [[partial derivative].sub.n][phi]. (29)

where [[partial derivative].sub.n][phi] is the boundary term introduced by discrete divergence acting on dual 2-cochain [*.sup.(1)] [[bar.d].sup.(0)][phi], and it is a dual 3-cochain in the whole space and a dual 2-cochain on the boundary surface.

A vectorial version of this, [nabla] x [nabla] x operator ([parallel]), can be written similarly. We will then examine how this operator applies to a vector field E(r). First, this vector field E(r) can be viewed as a primal 1-cochain E. Then the first curl operator acts as d (1) and will generate a primal 2-cochain:

[nabla] x E(r) [??] [[bar.d].sup.(1)] E. (30)

Then a hodge star operator is needed for the second curl operator to function appropriately, which leads to a dual 1-cochain [*.sup.(2)][[bar.d].sup.(1)]E. Then the second curl operator will introduce an extra boundary term as in Equation (20):

[nabla] x [nabla] x E(r) [??] [[([[bar.d].sup.(1)]).sup.T] [*.sup.(2)] [[bar.d].sup.(1)]] E + i[omega][mu][J.sub.s]. (31)

Here, coefficient i[omega][mu] is added to be consistent with Maxwell's equations, and Js represents the tangential component of [nabla] x E(r) or H(r) on the boundary.

3. WEDGE PRODUCT

In exterior calculus, wedge product is used to construct higher degree forms from lower degree ones. For example, with two 1-forms e and h in 3-D space, their wedge product e [conjunction] h is a 2-form. And this corresponds to the fact that the cross product of two vector fields E(r) and H(r) gives rise to a new vector field S(r) = E x H. Similarly, a wedge product between a 1-form e and a 2-form d in 3-D space gives rise to a 3-form e [conjunction] d, which refers to a scalar field generated by dot product between two vector fields. In fact, the dot product between electric field E(r) and displacement field D(r) is proportional to electric energy density, which can be integrated over a volume to obtain stored electric energy.

In this paper, we only need to introduce a discrete treatment of wedge product between a primal cochain and its corresponding dual cochain, which was described briefly in . With a primal fc-cochain b and a dual (n - k)-cochain h, their wedge product should give rise to a n- cochain, or a volume cochain. In practice, we are more interested in the integration of a volume form or a volume cochain, so it is not crucial whether this n-cochain should be a primal one or dual one.

Since cochains are column vectors, we use an inner product between two column vectors to represent the volume integration of their wedge product. For example, with a primal 1- cochain E representing the electric field and dual 2-cochain D for the electric displacement field in a given mesh, the stored electric energy (up to a constant difference) can be represented as:

[integral] [E(r) x D(r)] dV [??] [E.sup.T] x D (32)

Here, the dot product is valid since a dual 2-cochain have the same length with a primal 1-cochain, which is [N.sub.1], the number of primal edges.

Similarly, we can also represent the magnetic energy in a discrete fashion, which is an inner product between a primal 2-cochain B and a dual 1-cochain H both with length [N.sub.2]

[infinity] [B(r) x H(r)] dV [??] [B.sup.T] x H (33)

Moreover, we can also use this to represent the integration of two scalar functions, [phi](r) and [rho](r). We think [phi](r) and [rho](r) as a primal 0-cochain [phi] and dual 3-cochain [rho]. Then the integration of their product over a volume will be:

[integral] [[phi](r) [rho](r)] dV [??] [[phi].sup.T] x [rho]. (34)

Next we would like to show this treatment of discrete wedge product is rigorous under local constant field assumption on arbitrary 2-D surfaces. Since the relation in Equation (34) is trivial, we will focus on the wedge product between two vector fields. To show this, a mathematical identity needs to be introduced.

3.1. A Mathematical I d entity

In Figure 6, we use [mathematical expression not reproducible] to represent three unit vectors along three primal edges, and [l.sub.i], [s.sub.i] to denote the length of each primal and dual edge for simplicity. Then, below relation holds true:

[mathematical expression not reproducible] (35)

where S is the area of this triangle, and [bar.I] is a 2 x 2 identity matrix. We will prove this in the appendix.

Then with two local constant two-dimensional vector field E(r) and D(r), the integration of their dot product over the triangle in Figure 6 is

[mathematical expression not reproducible] (36)

where [E.sub.i] = (E x [[??].sub.i])[l.sub.i] is the value of cochain E on primal edge [l.sub.i], and [D.sub.i] = (D x [[??].sub.i])[s.sub.i] is the value of cochain D on dual edge [s.sub.i]. In a 2-D mesh with multiple triangles, the value of dual cochain D on each dual edge normally comes from two triangles, which can be combined. Then the relation becomes

[mathematical expression not reproducible] (37)

where [A.sub.i] refers to i-th primal face. Therefore, Equation (37) shows that the dot product between a primal 1-cochain and a dual 1-cochain is a good approximation of the integration of E x D. The relation in Equation (35) is also true with [[??].sub.i] changing to [mathematical expression not reproducible]

[mathematical expression not reproducible] (38)

With this relation, we can prove

[integral][B(r) x H(r)] dS [approximately equal to] [[N.sub.1].summation over (j=1)] (39)

in a 2-D mesh in the same way.

Moreover, we can also consider the integration of E x H on a 2-D surface, which is not necessarily a plane, as a wedge product

[mathematical expression not reproducible] (40)

where [??] is the normal unit vector of this surface.

It should be pointed out that in a general 3-D mesh, although the counterpart relation of Equation (35) does not hold exactly, we still adopt inner product in Equations (37) and (39) as the volume integration of discrete primal-dual wedge product. With this discrete treatment of wedge product, we can derive the discrete version of Huygens' principle, Poynting's theorem, and reciprocity theorem.

4. DISCRETE ELECTROMAGNETIC THEORY

With the discrete calculus introduced in Section 2, in 3-D space without a boundary, discrete Maxwell's equations can be written as

[[bar.d].sup.(1)] E = i[omega]B, (41a)

[([[bar.d].sup.(1)]).sup.T] H = -i[omega]D + J, (41b)

[[bar.d].sup.(2)] B = 0, (41c)

-[([[bar.d].sup.(0)]).sup.T] D = [rho], (41d)

where relations in Equation (15) have been applied.

For any k-cochain [[omega].sub.k], it is easy to verify that:

[[gamma].sub.k+2] = [[bar.d].sup.(k+1)][[bar.d].sup.(k)][[omega].sub.k] = 0, k [less than or equal to] 1 in a 3-D mesh, and k = 0 in a 2-D mesh.

In fact, this is because, in the boundary of the boundary of a k-simplex, a (k - 2)-simplex always appears twice with different signs . In exterior calculus language, this means that any exact form is automatically closed. The opposite way is also guaranteed by Poincare lemma. Therefore, in a 3-D mesh, due to relations [[bar.d].sup.(1)] x [[bar.d].sup.(0)] = 0 and [[bar.d].sup.(2)] x [[bar.d].sup.(1)] = 0, [nabla] x [nabla] = 0 and [nabla] x [nabla] x = 0 are exactly preserved. In fact, the second relation gives rise to charge continuity equation, also derived in .

[([[bar.d].sup.(0)]).sup.T] x [([[bar.d].sup.(1)]).sup.T] H = 0 [??] [([[bar.d].sup.(0)]).sup.T] (-i[omega]D + J) = 0 [??] i[omega]p + [([[bar.d].sup.(0)]).sup.T] J = 0 (42)

This means that Maxwell's equations in DEC are consistent with charge conservation, which is the same with continuous theory. Therefore, numerical simulations based on DEC will not generate spurious charges.

Moreover, two constitutive relations are needed to complete Maxwell's equations

[mathematical expression not reproducible]. (43)

Here, two Hodge star operators [mathematical expression not reproducible] map primal cochains E and B to dual cochains D and H, superscripts refer to the degree of cochains they are acting on, and subscripts e and fj.-1 means that they are characterized by these two material information, which is described in detail in . For homogeneous medium, [mathematical expression not reproducible].

Noted that the choice of E and B as primal cochains and H, D and J as dual cochains is only by convention, not a requirement. We have shown that the other choice is beneficial for PEC boundary condition implementation .

4.1. Green's Function

4.1.1. Scalar Waves

In continuous theory, a Green's function of a wave equation is the solution of the wave equation with a point source. When this point source solution is known, we can obtain the solution to a general source through superposition. For example, the scalar wave equation in 3-D can be written as:

([[nabla].sup.2] + [k.sup.2])[psi](r) = s(r). (44)

Then by treating [psi](r), s(r) as a primal 0-cochain [psi] and a primal 3- cochain s (both are [N.sub.0] x 1 column vectors), we can rewrite this equation in a discrete fashion as:

[([[bar.d].sup.(0)]).sup.T] [*.sup.(1)][[bar.d].sup.(0)] - [k.sup.2][*.sup.(0)]] x [psi] = -s. (45)

Noted that the space considered is assumed to be infinite, so that the boundary term is dropped. The scalar Green's function g(r, r') is the solution to following equation:

([[nabla].sup.2] + [k.sup.2])g(r, r') = -[delta](r - r'). (46)

Here, for each given r', g(r, r') should be treated as a [N.sub.0] x 1 column vector (0-cochain). Then with [N.sub.0] possible r', Green's function g(r, r') should be considered as a [N.sub.0] x [N.sub.0] matrix g. Similarly, for any given r', [delta](r - r') should be considered as a [N.sub.0] x 1 column vector (3-cochain) with only one nonzero entry with value 1, where r = r'. Therefore, a discrete [delta](r - r') can be replaced by a [N.sub.0] x [N.sub.0] identity matrix.

[[([[bar.d].sup.(0)]).sup.T] [*.sup.(1)][[bar.d].sup.(0)] - [k.sup.2][*.sup.(0)] x [bar.g] = [bar.I]. (47)

Here the boundary term is also dropped. Then,

[bar.g] = -[[[([[bar.d].sup.(0)]).sup.T] [*.sup.(1)] [[bar.d].sup.(0)] + [k.sup.2][*.sup.(0)]].sup.-1] (48)

is the discrete solution of Green's function in an infinite homogeneous space. Notice that since Hodge operators are diagonal, this discrete Green's function is a symmetric matrix. This coincides with the fact that for the continuous solution, g(r, r') = g(r', r).

Then we left multiply Equation (45) with [[bar.g].sup.T], and right multiply Equation (47)'s transpose with [psi]:

[[bar.g].sup.T] x (45) [??] x [[([[[bar.d].sup.(0)]).sup.T] [*.sup.(1)] [[bar.d].sup.(0)] - [k.sup.2] [*.sup.(0)]] x [psi] = -[[bar.g].sup.T] x s (49a)

[(47).sup.T] x [psi] [??] [[bar.g].sup.T] x [[([[bar.d].sup.(0)]).sup.T] [*.sup.(1)] [[bar.d].sup.(0)] - [k.sup.2] [*.sup.(0)]] x [psi] = [psi] (49b)

Therefore, we can represent the solution [psi] due to an arbitrary source s with Green's function as:

[psi] = -[[bar.g].sup.T] x s. (50)

Obviously, this is a discrete analog of the continuous version

[psi] (r) = -[[integral].sub.V] dr'g (r, r') s (r'). (51)

In fact, this can also be understood from wedge product between a primal 0- cochain and a dual 3-cochain. More specifically, for any given r, scalar field g(r, r'), with respect to r', corresponds to a primal 0-cochain, and s(r') is viewed as a dual 3-cochain. In Section 3, we have introduced the volume integration of a wedge product between these two cochains can be represented by an inner product. Therefore, with [N.sub.0] possible r, this integration can be represented by a matrix vector product.

4.1.2. Vector Waves

For vector waves, a dyadic Green's function is introduced to describe the solution to a Hertzian dipole source. First, the equation for the vector wave in a homogeneous, isotropic medium is

[nabla] x [nabla] x E(r) - [k.sup.2]E(r) = i[omega][mu]J(r). (52)

where k = [omega][square root of ([mu][epsilon])] is the wave number.

Then, similarly, by treating E(r) and J(r) as a primal 1-cochain E and a dual 2- cochain J (both are [N.sub.1] x 1 column vectors), this vector wave equation can be written in a discrete fashion as:

[[([[bar.d].sup.(1)]).sup.T] [*.sup.(2)] [[bar.d].sup.(1)] - [k.sup.2][*.sup.(1)] ] x E = i[omega][mu]J. (53)

Here, the equation is assumed for the whole space, so boundary term is also dropped.

The dyadic Green's function [bar.G] (r, r') is solution to following equation:

[nabla] x [nabla] x [bar.G] (r, r') - [k.sup.2] [bar.G] (r, r') = [bar.I] [delta] (r - r'). (54)

Here, slightly different from the scalar case, for each given r', [bar.G](r, r') is a tensor field, which is a 3 x 3 matrix at any r. In fact, for each given Hertzian dipole current source placed at primal edge around r', if the unit vector along this primal edge is [??]', then [bar.G] (r, r') x [??]' is a vector field and can be treated as a primal 1-cochain, which is a [N.sub.1] x 1 ([N.sub.1] is the number of primal edges) column vector. Then with [N.sub.1] possible r', [bar.G] (r, r') can be viewed as an [N.sub.1] x [N.sub.1] matrix, and [bar.I][delta] (r - r') on the right hand side should be treated as an [N.sub.1] x [N.sub.1] identity matrix. Then, with DEC, the equation for dyadic Green's function can be written as:

[[([[bar.d].sup.(1)] [*.sup.(2)] [[bar.d].sup.(1)] - [k.sup.2][*.sup.(1)]).sup.T] ] x [bar.G] = [bar.I]. (55)

Therefore,

[bar.G] = [[[([[bar.d].sup.(1)]).sup.T] [*.sup.(2)] [[bar.d].sup.(1)] - [k.sup.2][*.sup.(1)]].sup.-1] (56)

is the discrete dyadic Green's function in homogeneous, isotropic medium. Because Hodge star matrices are diagonal, hence symmetric, this [bar.G] is also a symmetric matrix. Combining Equations (53) and (55), we can obtain

E = i[omega][mu][[bar.G].sup.T] x J. (57)

Clearly, this is the discrete version of

E(r) = i[omega][mu] [[integral].sub.V] [bar.G] (r, r') x J (r') dr'. (58)

We can also view Equation (57) as a discrete wedge product. To see this, we can write the value of primal 1-cochain E on a primal edge [l.sub.i] around r by using Equation (58):

[mathematical expression not reproducible] (59)

For a given primal edge [mathematical expression not reproducible] gives rise to a vector field, and can be viewed as a primal 1-cochain. In fact, this primal 1-cochain is just the i-th column of [bar.G]. Therefore, the integration in Equation (58) can be viewed as a wedge product between this primal 1-cochain and dual 2-cochain J. From the discussion in Section 3, the integration of this wedge product can be carried out through a vector inner product as in Equation (57).

So far, we have introduced the Green's function for both scalar waves and vector waves. Although we assumed homogeneous medium for simplicity, inhomogeneous medium can be considered similarly with characterized Hodge star operators, which we discussed in . Next we will include the boundary effects and obtain Huygens' principle for both scalar waves and vector waves.

4.2. Huygens' Principle

Huygens' principle states that the wave field on a closed surface S can determine the field elsewhere. This principle is valid for both scalar and vector waves. We will derive the discrete version of this principle in the following. The scalar wave case is considered first.

4.2.1. Scalar Waves

When deriving Equations (45) and (47), we dropped the boundary term since the whole space is considered. However, generally, the interested space is a finite region. Hence the boundary term cannot be ignored.

As shown in Figure 7, the considered region is bounded by surface S and [S.sub.inf], and the source only exists within S. Then, by using the conclusion in Equation (29), the discrete scalar wave equation with a boundary term can be written as:

[[([[bar.d].sup.(0)]).sup.T] [*.sup.(1)] [[bar.d].sup.(0)] - [k.sup.2][*.sup.(0)]] x [psi] - [[partial derivative].sub.n][psi] = 0. (60)

Here, since [S.sub.inf] is located at infinity, its contribution is assumed to be zero. So [[partial derivative].sub.n][psi] has nonzero value only on S. More specifically, [[partial derivative].sub.n][psi] is a dual 3- cochain with nonzero values only on the dual 3-cells associated with vertices on S. In fact, as mentioned, [[partial derivative].sub.n][psi] can also be viewed as a dual 2- cochain on S, which is a 2-D manifold. The value of [[partial derivative].sub.n][psi] on the dual 2-cell [mathematical expression not reproducible] on S in Figure 7 is

[mathematical expression not reproducible]. (61)

Here, the last approximation is given because this form is easier to implement boundary conditions. The equation for scalar Green's function is written as:

[[([[bar.d].sup.(0)]).sup.T] [*.sup.(1)] [[bar.d].sup.(0)] - [k.sup.2] [*.sup.(0)]] x [bar.g] - [[partial derivative].sub.n][bar.g] = [bar.I]. (62)

Here, [[partial derivative].sub.n][bar.g] is a [N.sub.0] x [N.sub.0] matrix denoting the boundary term, and its (i, j) element is defined as:

[mathematical expression not reproducible] (63)

where [r.sub.j] is the j-th primal vertex, and [A.sup.[partial derivative].sub.i] is dual face on the boundary associated with i-th primal vertex [r.sub.i], such as [mathematical expression not reproducible] in Figure 7. With [N.sup.[partial derivative].sub.0] vertices on S, [[partial derivative].sub.n][bar.g] will only have [N.sup.[partial derivative].sub.0] nonzero rows. It is noted that each column of [[partial derivative].sub.n][bar.g] can be viewed as a dual 3-cochain, or a dual 2-cochain on S, just as [[partial derivative].sub.n][psi].

Combining Equations (60) and (62), we can derive

[psi] = [[bar.g].sup.T] x [[partial derivative].sub.n][psi] - [([[partial derivative].sub.n][bar.g]).sup.T] x [psi] (64)

Here, the two matrix-vector production can be understood as a series of wedge product between primal 0-cochains and dual 3-cochains as defined in Equation (34), or a wedge product between primal 0-cochains and dual 2-cochains on 2-D manifold S, which refers to a surface integral. Clearly, this is a discrete analog to the continuous theory:

[mathematical expression not reproducible]. (65)

Moreover, if we set g (r, r') = 0 for r [member of] S, then the first term on the right hand side of Equation (64) will vanish and become

[psi] = -[([[partial derivative].sub.n][bar.d]).sup.T] x [psi]. (66)

On the other hand, if we require [??] x [nabla] g (r, r') = 0 for r [member of] S, then instead, the second term of (64) will disappear and

[psi] = [[bar.g].sup.T] x [[partial derivative].sub.n][psi]. (67)

Equations (64), (66) and (67) are various forms of discrete Huygens' principle for scalar field depending on the definitions of g(r, r'). This principle states that only the value of [psi] or [[partial derivative].sub.n][psi] on surface S are needed to obtain the field value at any point between S and [S.sub.inf].

4.2.2. Electromagnetic Waves

For electromagnetic waves, we consider the same geometry, as shown in Figure 8. It is noted that the red solid line [mathematical expression not reproducible], connecting [d.sub.0] and [d.sub.1], is from [d.sub.0] to middle point of primal edge [bar.[p.sub.0][p.sub.1]] then to [d.sub.1], and not necessarily straight. And dual face [mathematical expression not reproducible] is perpendicular to [bar.[p.sub.0][p.sub.1]], and is inside the mesh.

Including a boundary term induced from the Laplacian operator, shown in Equation (31) in Section 2, in a source-free region bounded by surface S and [S.sub.inf], Equation (53) is adjusted as

[[([[bar.d].sup.(1)]).sup.T] [*.sup.(2)][[bar.d].sup.(1)] - [k.sup.2][*.sup.(1)]] x E + i[omega][mu][J.sub.s] = 0. (68)

The contribution of [S.sub.inf] is also assumed to be zero, then [J.sub.s] has nonzero value only on S. More specifically, [J.sub.s] is a dual 2-cochain with nonzero entities only on the incomplete dual faces at the boundary, e.g., [mathematical expression not reproducible] in Figure 8. Due to the one to one relation between dual faces and primal edges, with [N.sup.[partial derivative].sub.1] primal edges on S, [J.sub.s] would have [N.sup.[partial derivative].sub.1] nonzero entries. For example, the value of [J.sub.s] defined on the dual face [mathematical expression not reproducible] is

[mathematical expression not reproducible] (69)

It should be noted that [J.sub.s] is a dual 1-cochain on 2-D manifold S, since it is defined on dual edges on S such as [mathematical expression not reproducible]. By writing as the second line of Equation (69), we can easily construct the wedge product between [J.sub.s] and a primal 1-cochain on S, as we derived in Equation (40) in Section 3.

The third line of Equation (69) is listed because this form is easier to implement certain kind of boundary condition, such as absorption boundary conditions (ABCs) and impedance boundary condition (IBC). For example, the first-order ABC reads as

[mathematical expression not reproducible].

Then, the line integral of [??] x H in Equation (69) can be replace by the value of primal cochain E on [bar.[p.sub.0][p.sub.1]] .

Similarly, with a boundary term, the equation for vector Green's function becomes

[[([[bar.d].sup.(1)]).sup.T] [*.sup.(2)] [[bar.d].sup.(1)] - [k.sup.2][*.sup.(1)]] x [bar.G] + [bar.K] = [bar.I]. (70)

Here, [bar.K] is an [N.sub.1] x [N.sub.1] matrix with only [N.sub.1,[partial derivative]] nonzero rows, similar with [[partial derivative].sub.n][bar.g] introduced for scalar wave. The (i,j) element [[[bar.K]].sub.i,j] associated with j-th primal edge and i-th dual face at surface, which only has one edge [L.sup.[partial derivative].sub.i] on the surface, is defined as:

[mathematical expression not reproducible] (71)

It is noted that with the j-th primal edge [mathematical expression not reproducible] gives rise to a vector field such as H(r). The tangential component of this vector field on the boundary, as defined in Equation (71), should be viewed as a dual 2-cochain in a 3-D mesh, or a dual 1-cochain on S, which is an [N.sub.1] x 1 column vector with [N.sup.[partial derivative].sub.1] nonzero entries, just as [J.sub.s]. Therefore, with [N.sub.1] possible [l.sub.j], these cochains form an [N.sub.1] x [N.sub.1] matrix, which is [bar.K].

Combining Equations (68) and (70), we can obtain the vectorial Huygens' principle:

E = [[bar.K].sup.T] x E - i[omega][mu][[bar.G].sup.T] x [H.sub.[parallel]] (72)

This formula can also be understood as a discrete wedge product between dual 1- cochains and primal 1-cochains on surface S. More specifically, E and each column of [bar.G] are primal 1-cochains, while [J.sub.s] and each column of [[bar.G].sub.[parallel]] are dual 1-cochains on S. Therefore, an inner product between them refers to a surface integral. Equation (72) is a discrete analog of the continuous theory 

[mathematical expression not reproducible]. (73)

The extra negative sign is due to the definition of normal unit vector [??].

4.3. Poynting's Theorem

The complex Poynting's theorem is written as 

[nabla] x (E x [H.sup.*]) = i[omega][H.sup.*] x B - i[omega]E x [D.sup.*] - E x [J.sup.*], (74a)

[mathematical expression not reproducible]. (74b)

where Equation (74b) is the integral form of Poynting's theorem on domain [OMEGA] and its boundary [partial derivative][OMEGA].

In Section 3, we have introduced the way to represent the right hand side of Equation (74b) using discrete wedge product as

RHS of (74b) [??] i[omega] [([H.sup.*]).sup.T] x B - i[omega][E.sup.T] x [D.sup.*] - [E.sup.T] x [J.sup.*] (75)

Invoking the discrete Ampere's law in Equation (41b) with a boundary,

-i[omega]D + J = [([[bar.d].sup.(1)]).sup.T] H + [J.sub.s]. (76)

Equation (75) can be adjusted as

i[omega] [([H.sup.*]).sup.T] x B - [E.sup.T] x [J.sup.*] - [E.sup.T] x [[([[bar.d].sup.(1)]).sup.T] [H.sup.*] + [J.sup.*.sub.s] - [J.sup.*]] = -[E.sup.T] x [J.sup.*.sub.s] (77)

Here we have applied [[bar.d].sup.(1)] E = i[omega]B. Next we need to show that the left hand side of Equation (74b) is represented by -[E.sup.T] x [J.sup.*.sub.s]. In fact, E and [J.sub.s] are primal 1-cochain and dual 1-cochain on two- dimensional manifold [partial derivative][OMEGA], representing the tangential components of E(r) and H(r) on the surface. In other words, [J.sub.s] can also be written as a dual 1-cochain [H.sub.[parallel]] defined on the 2-D surface. For example, the value of [H.sub.[parallel]] on [mathematical expression not reproducible] in Figure 8 is defined as

[mathematical expression not reproducible] (78)

Then the discrete counterpart of the left hand side of Equation (74b) can be represented by the inner product of E and [J.sub.s] or [H.sub.[parallel]], as shown in Equation (40). The negative sign is due to the contradiction between the direction of dual edges on boundary, and normal direction of associated dual faces. For example, in Figure 8, if the direction of [bar.[p.sub.0][p.sub.1]] is assumed as [??], which is from [p.sub.0] to [p.sub.1]. Then, the normal direction of dual face [mathematical expression not reproducible] is also [??], while the direction of dual edge [mathematical expression not reproducible] on the surface is [mathematical expression not reproducible], which is from [d.sub.1] to [d.sub.0] and contradicts with the normal direction of [mathematical expression not reproducible].

Therefore, we have shown that the discrete Poynting's theorem still holds true and is written as

[mathematical expression not reproducible] (79)

Here, the constitutive relations introduced in Equation (43) are applied, and the superscripts of Hodge star operators are dropped for simplicity. With a lossless and source free region, Hodge star operators are real, then

[mathematical expression not reproducible]. (80)

This means that no net time-average power-flow out of or into domain O.

4.4. Reciprocity Theorem

A medium is called reciprocal medium when its permittivity [bar.[epsilon]] and permeability [bar.[mu]] are symmetric tensors. In a reciprocal medium, with an electric field [E.sub.1] generated by [J.sub.1] and [E.sub.2] generated by [J.sub.2], they should satisfy

<[E.sub.2], [J.sub.1]> = <[E.sub.1], [J.sub.2]>. (81)

We will prove this theorem within DEC only assuming two Hodge star operators [*.sup.(1).sub.[epsilon]] and [mathematical expression not reproducible] are symmetric. In fact when the medium is isotropic, these two Hodge star operators are diagonal with DEC. With anisotropic medium, they are still symmetric, but not diagonal any more, which will be shown in our future work.

Since E field and source J can be represented by a primal 1-cochain and dual 2- cochain, respectively, the inner product in Equation (81) can be viewed as a wedge product between cochains. In other words, we need to prove

[E.sup.T.sub.2] x [J.sub.1] = [E.sup.T.sub.1] x [J.sub.2]. (82)

Cochains [E.sub.1] and [E.sub.2] are generated by source cochains [J.sub.1] and [J.sub.2] via equations

[mathematical expression not reproducible], (83a)

[mathematical expression not reproducible]. (83b)

Therefore,

[mathematical expression not reproducible]. (84)

Noted that this relationship is exact in this discrete space.

5. CONCLUSIONS

We have shown that by using DEC and a discrete wedge product, discrete electromagnetic theory can be built on an arbitrary simplicial mesh. The fact that theorems of electromagnetic theory are still preserved in this discrete world, shows that discrete electromagnetic theory with DEC is self-consistent and self-contained. Therefore, these discrete theorems can be used to check the correctness of numerical simulations. Moreover, the discrete Green's functions can be used to solve complex scattering problems, which will be explored in our future works.

ACKNOWLEDGMENT

We would like to acknowledge the following funding sources: UIUC CAS, AF Sub RRI P00539, NSF ECCS 169195, ANSYS Inc P037497 and the George and Ann Fisher Professorship.

APPENDIX A. PROOF OF EQUATION (35)

Since this is a 2 x 2 matrix equation, we only need to prove that both sides of Equation (35) acting on an arbitrary 2-D vector give rise to the same result. Since any vector on a 2-D plane can be written as a linear supposition [of [??].sub.1] and [[??].sub.2], we can only use these two vectors to prove our point.

RHS of (35) x [[??].sub.1] = S[[??].sub.1] (A1)

[mathematical expression not reproducible] (A2)

For any triangle, the following relation holds true

[mathematical expression not reproducible] (A3)

Therefore,

[mathematical expression not reproducible] (A4)

[mathematical expression not reproducible] (A5)

To prove that the right hand side of Equations (A2) and (A5) are identical, we need to prove the following two relations

S = [l.sub.1] [[s.sub.1] - [s.sub.3] ([[??].sup.T.sub.3] - [[??].sub.1])] (A6)

[mathematical expression not reproducible] (A7)

As shown in Figure A1, we have the following relations

[[??].sup.T.sub.2] x [[??].sub.1] = -cos ([[theta].sub.2]) (A8)

[[??].sup.T.sub.3] x [[??].sub.1] = -cos ([[theta].sub.1]) (A9)

[mathematical expression not reproducible] (A10)

[mathematical expression not reproducible] (A11)

Therefore, Equations (A6) and (A7) are proved, which leads to

LHS of (35) x [[??].sub.1] = RHS of (35) x [[??].sub.1] (A12)

Since the same derivation also works for [[??].sub.2],

LHS of (35) x [[??].sub.2] = RHS of (35) x [[??].sub.2] (A13)

So Equation (35) is proved.

Caption: Figure A1. Geometry of an arbitrary triangle. Here, [d.sub.0] is the circumcenter, [l.sub.1], [l.sub.2], [l.sub.3] are three primal edges, and [d.sub.1], [d.sub.2], [d.sub.3] are three dual edges.

REFERENCES

[1.] Chew, W. C., Waves and Fields in Inhomogeneous Media, IEEE Press, New York, 1995.

[2.] Jin, J.-M., The Finite Element Method in Electromagnetics, John Wiley & Sons, 2015.

[3.] Chew, W. C., "Electromagnetic theory on a lattice," Journal of Applied Physics, Vol. 75, 4843-4850, 1994.

[4.] Deschamps, G. A. "Electromagnetics and differential forms," Proceedings of the IEEE, Vol. 69, No. 6, 676-696, 1981.

[5.] Bossavit, A., "Differential forms and the computation of fields and forces in electromagnetism," Eur. J. Mech. B, Vol. 10, No. 5, 474-488, 1991.

[6.] Warnick, K. F., R. H. Selfridge, and D. V. Arnold, "Electromagnetic boundary conditions and differential forms," IEE Proceedings Microwaves Antennas and Propagation, Vol. 142, No. 4, 326-332, 1995.

[7.] Warnick, K. and D. Arnold, "Green forms for anisotropic, inhomogeneous media," Journal of Electromagnetic Waves and Applications, Vol. 11, No. 8, 1145-1164, 1997.

[8.] Teixeira, F. L. and W. C. Chew. "Lattice electromagnetic theory from a topological viewpoint," Journal of Mathematical Physics, Vol. 40, No. 1, 169-187, 1999.

[9.] Teixeira, F. L., "Differential forms in lattice field theories: An overview," ISRN Mathematical Physics, Vol. 2013, 2013.

[10.] Desbrun, M., A. N. Hirani, and M. Leok, "Discrete exterior calculus," arXiv preprint math/0508341, 2005.

[11.] Hirani, A. N., "Discrete exterior calculus," Ph.D. Thesis, California Institute of Technology, 2003.

[12.] Desbrun, M., K. Eva, and Y. Tong, "Discrete differential forms for computational modeling," Discrete Differential Geometry, 287-324, Birkhauser Basel, 2008.

[13.] Hirani, A. N., K. Kalyanaraman, and E. B. VanderZee, "Delaunay hodge star," Computer-Aided Design Vol. 45, No. 2, 540-544, 2013.

[14.] Chen, S. and W. C. Chew, "Numerical electromagnetic frequency domain analysis with discrete exterior calculus," arXiv preprint arXiv:1704.05145, 2017.

[15.] Madsen, N. K. and R. W. Ziolkowski, "A three-dimensional modified finite volume technique for maxwell's equations," Electromagnetics, Vol. 10, 147-161, 1990.

[16.] Clements, M. and T. Weiland, "Discrete electromagnetism with the finite integration technique," Progress In Electromagnetics Research, Vol. 32, 65-87, 2001.

[17.] Na, D.-Y., H. Moon, Y. A. Omelchenko, and F. L. Teixeira, "Local, explicit, and charge-conserving electromagnetic particle-in-cell algorithm on unstructured grids," IEEE Transactions on Plasma Science, Vol. 44, No. 8, 1353-1362, 2016.

[18.] Chen, S. and W. C. Chew, "Discrete electromagnetic theory with exterior calculus," PIERS Proceedings, 896-897, Shanghai, China, August 8-11, 2016.

Shu C. Chen (1) and Weng C. Chew (2), *

Received 15 May 2017, Accepted 3 July 2017, Scheduled 9 July 2017

* Corresponding author: Weng Cho Chew (w-chew@illinois.edu).

(1) Department of Physics, University of Illinois at Urbana-Champaign, USA. 2 Department of Electrical and Computer Engineering, University of Illinois at Urbana-Champaign, USA.

([dagger]) Generally, the electric displacement field D(r) is treated as a dual cochain. Here we only use it to illustrate the idea since the divergence of magnetic field B(r) is always zero. ([double dagger]) Usually, dual faces or volumes are not simplices, as in Figure 3. So we use dual cells to denote them.

([section]) In this paper, we use A to denote dual faces, which are usually polygons.

([double dagger]) Usually, dual faces or volumes are not simplices, as in Figure 3. So we use dual cells to denote them.

([parallel]) Strictly speaking, this is not the Laplacian operator for a vector field, -W- needs to be added.

Caption: Figure 1. A primal face [S.sub.i] with oriented edges. Here [[??].sub.i] is assumed as the positive direction of [S.sub.i].

Caption: Figure 2. (a) Primal face [S.sub.i] in a 2-D mesh. (b) Primal volume [V.sub.i] in a 3-D mesh.

Caption: Figure 3. Dual cells at the boundary (bold edges) in a 2-D mesh. Here, [d.sub.1], ..., [d.sub.6] are circumcenters, and [e.sub.1], [e.sub.2] are edge centers.

Caption: Figure 4. Dual cells at the boundary (bold edges) in 2-D space.

[??] direction is out of the paper.

Caption: Figure 5. Dual face at the boundary (surrounded by dashed edges) in 3-D space. [??] is a unit vector normal to the surface, and [mathematical expression not reproducible] (orange region) is the dual 2-cell on the surface boundary associated with [p.sub.0]. Caption: Figure 6. A triangle with primal edges [l.sub.1], [l.sub.2], [l.sub.3] and dual edges [s.sub.1], [s.sub.2], [s.sub.3]. Here, [d.sub.0] is the circumcenter.

Caption: Figure 7. The geometry for Huygens' principle derivation. [p.sub.0] is a vertex on the surface S, [??] is the unit normal vector. The orange shaded region [mathematical expression not reproducible] is obtained by connecting circumcenters of triangles around [p.sub.0] on S.

Caption: Figure 8. The geometry for vector Huygens' principle derivation. Here, [bar.[p.sub.0][p.sub.1]] is a primal edge on surface S. [mathematical expression not reproducible] is the dual face associated with [bar.[p.sub.0][p.sub.1]], [d.sub.0] and [d.sub.1] are the circumcenters of two triangles containing [bar.[p.sub.0][p.sub.1]], and [mathematical expression not reproducible] is the red solid line connecting [d.sub.0] and [d.sub.1]]