# A new element-oriented model for computational electromagnetics.

I. INTRODUCTION

To solve complex real-world electromagnetic problems, several methods have been developed over the past three decades taking advantage of the powerful computer resources [1-3]. Each one of them has advantages/disadvantages over the others. Some of the most common methods can be classified as follows:

* the finite difference method (FDM),

* the finite integral method (FIT),

* the finite element method (FEM),

* the finite volume method (FVM),

* the method of moments (MoM), and

* the boundary element method (BEM).

All these numerical methods (and their extensions) have been developed separately, but they have numerous commonalities [46]. Numerical methods yield approximate solutions to the governing equations through the discretization of space and time. If the user of a model is unaware of the details of the numerical method, including the derivative approximations, the scale of discretization, and the matrix solution techniques, significant errors can be introduced and remain undetected [8].

In this paper, we initiate the development of a new Four-dimensional numerical model completely oriented object that we call an Element-Oriented Model (EOM). Its general formulation is built on new conceptual views. Therefore, it can be used to review and to perform the most known methods. We consider the Element-Oriented Modeling in electromagnetism as a natural process to systematize the electromagnetic concepts in a model by distinguishing objects or individual elements that include information (state or data values) and functionality (behavior). Using an object-oriented approach to arranging a set of equations allows us to group particular physical entities together with common functionalities or actions associated with that information. The state of an object encompasses all of the (usually static) properties of the object plus the current (usually dynamic) values of each of these properties. Through the Object-Oriented Paradigm, we can see electromagnetic problems as a collection of interacting objects. In practice, this approach is based on the strategies and techniques (often called Object-Oriented analysis) for establishing a physical model that distinguishes objects, their properties and their actions from the perspective of the classes and objects founded in the vocabulary of the electromagnetism theory. In this paper, we investigate the differential forms and the concept of topological laws in electromagnetism to formulate a general Element-Oriented Model. Notice that topology was recognized by Leibniz, Gauss and Maxwell to play a central role in the electromagnetism theory [11]. This tool is not largely exploited in computational physics [10]. The advantages of the topological approach are various:

* It provides an important framework connection between the electromagnetic fields and the geometry [7-11,15,17,23,2529,33,35].

* It provides an excellent viewpoint for the discretization of electromagnetic equations [7, 9, 12-15, 35].

* It gives discrete formulations preserving many physically significant properties of the original problem [8-10,27,35].

* It permits to draw a dynamic representation of Maxwell's equations [9, 35].

* Furthermore, this approach facilitates the development of new numerical methods [9, 35].

Throughout this paper, we investigate this approach, and we propose an Element-Oriented Model (EOM) enjoying all the benefits of the Object-Oriented Programming (OOP) (abstraction, encapsulation, modularity, hierarchy, polymorphism and persistence). The first task is to find all objects with needed properties and behaviors. Then, we present a general finite formulation for this model. From which, we derive a Flexible Finite-Difference-Time-Domain (Flexible FDTD) algorithm. This formulation generalizes both the standard S-FDTD and the nonstandard NS-FDTD methods. To show its efficiency, we present a new nonstandard scheme derived directly from it. Through a numerical example, we show that this new scheme is more accurate than the S-FDTD and the known NS-FDTD methods. Because the S FDTD and the NS-FDTD methods present distinct discrete derivative operators, that they are chosen in this work. Indeed, many other methods like ADI-FDTD [18], LOD-FDTD [19], and CN FDTD [20] use the same discrete derivative operator as the S-FDTD method.

2. THE ELEMENT-ORIENTED MODEL

2.1. Geometric Objects

The concept of the geometric objects is ubiquitous in physical field theories [6-9, 17, 33]. In three-dimensional space, the integral form of Maxwell-Heaviside equations reveals four basic spatial geometrical objects such as points P, lines L, surfaces S and volumes V, and two temporal geometrical elements such as instants T and time-intervals I. However, this integral formulation does not clearly reflect the orientation of these objects neither does the axial/polar nature of electromagnetic fields associated with them [35]. Indeed, all geometric objects are endowed with two kinds of orientation: internal or external. We will conventionally add a tilde to distinguish externally oriented objects, from internally oriented ones, thus writing [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for externally oriented point, line, surface, volume, instant, and interval respectively. For more details on these geometric objects, we refer the reader to [8, 9, 33].

If we adopt a strict space-time viewpoint, we must consider space and time as one four-dimensional space. Therefore, the space-time objects can be considered as Cartesian products of a space object by a time object. So, we obtain a new space-time structure characterized by two cell complexes:

* A primal cell Complex K defined by:

K = {P, L, S, V, H} = {P, L, S, V} x {T, I} (1)

* A dual cell Complex [??] defined by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

The space-time objects in the primal cell complex K have inner orientation, whereas the space-time objects in the dual cell complex [??] have outer orientation. To distinguish space-time objects from merely spatial ones, we will use the symbols P, L, S, V, and H for the former and the symbols P, L, S, and V for the latter. H = V x I is an hypervolume in the cell complex K. On these geometric objects, we can apply the boundary operator d. So, we get, as examples [partial derivative]S = L; [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

These geometric objects are considered basic tools of algebraic topology. They are called oriented p-dimensional cells (p = 0, 1, 2, 3, 4) in a four-dimensional space (space-time as example). A generic sub domain of p-dimensional cells is called a p-chain. So, the boundary operator [partial derivative] can be shown as a linear mapping of the space of p-chains into that of (p - 1)-chains.

2.2. Physical Global Quantities of Electrodynamics

After defining the oriented geometric objects, the next task is to find their association with electromagnetic entities. For this, we define the following physical variables on the space-time structure characterized by K and [??]:

* [[PHI].sup.2] = ([[PHI].sup.e], [[PHI].sup.b]) is the electromagnetic flux associated to primal surfaces S. [[PHI].sup.2] is a pairing of [[PHI].sup.e] with [[PHI].sup.b] where [[PHI].sup.e] is the electric flux associated with the object [partial derivative]S x I and [[PHI].sup.b] is the magnetic flux associated with S x [partial derivative]I.

* [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the electromagnetic charge associated with dual volumes [??]. [[??].sup.3] is a pairing of [[??].sup.[rho]] with [[??].sup.j] where [Q.sup.[rho]] is the electric charge content associated with the object [??] x [??] and [Q.sup.j] is the electric charge flow associated with [??] x [??].

* [U.sup.1] = ([U.sup.a], [U.sup.v]) is the electromagnetic potential associated with the object L. [U.sup.1] is a pairing of [U.sup.a] with [U.sup.v] where [U.sup.a] is the potential associated with the object L x [partial derivative]I and [U.sup.v] is the electric potential associated with the object [partial derivative]L x I.

* [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the electromagnetic charge-current potential (as called in [27]) associated with the object [??]. [[??].sup.2] is a pairing of [[??].sup.d] with [[??].sup.h] where [[??].sup.d] is the electric flux associated with the object [??] x [partial derivative][??] and [[??].sup.h] is the magnetic voltage impulse (as called in [17]) associated with d[??] x [??].

Each Physical Global Quantity is indexed by a prime to emphasize the dimension of the object associated with it. Association of these quantities with geometric objects makes more understandable any electromagnetic model.

These global quantities are also called global variables by other authors such as in [17]. By integration of field functions (E, B, H, D, [rho], J) on space domains (lines, surfaces, volumes) and on time intervals. Some authors deduce all these global variables [8, 9, 33]. In the next section, we use a formal analysis based on the differential forms in electromagnetism and some basic tools of algebraic topology. We deduce global quantities and we establish the discrete topological equations.

2.3. Discrete Topological Equations

In terms of differential forms, Maxwell equations can be written as [30, 31]:

dE = - [partial derivative]B/[partial derivative]t; dD = [rho] (3)

dH = [partial derivative]D/[partial derivative]t + J; dB = 0 (4)

where

* d denotes the spatial exterior derivative; when applied to 0, 1 or 2-forms respectively, d is equivalent to grad, curl and div operators of vector calculus. When acting on a k-form, d produces a (k + 1)-form.

* E is the electric intensity 1-form.

* B is the magnetic flux density 2-form.

* J is the current density 2-form.

* [rho] is the charge density 3-form.

These equations are independent from the space metric. So, they are invariant under diffeomorphisms. The metric equations also called constitutive equations can be defined as [8,9,13,14,17,22,25]:

D = [*.sub.[epsilon]]E; B = [*.sub.[mu]]H; J = [*.sub.[sigma]]E (5)

where * is the spatial Hodge star operator.

In the four-dimensional space-time, we define some electromagnetic entities. Each entity is indexed by a prime to emphasize the dimension of its form. So, we define the electromagnetic strength F 2-form as:

[F.sup.2] = [E.sup.1] [conjunction] dt + [B.sup.2] (6)

Its dual * [F.sup.2] is expressed as:

*[F.sup.2] = ([*.sub.[mu]][B.sup.2]) [conjunction] dt [*.sub.[epsilon]][E.sup.1] (7)

*[F.sub.2] = [H.sup.1] [conjunction] dt - [D.sup.2] (8)

where [conjunction] is the exterior product and * is the four-dimensional Hodge star operator.

We define the electromagnetic charge J as:

[J.sup.3] = [[rho].sup.3] - [J.sup.2] [conjunction] dt (9)

We define the electromagnetic potential [epsilon] as:

[[epsilon].sup.1] = [V.sup.0] [conjunction] dt + [A.sup.1] (10)

where V is the electric potential 0-form and A is the vector potential 1-form.

By applying the four-dimensional exterior derivative d to the Equations (6), (8), (9) and (10), we get:

d[F.sup.2] = [0.sup.3] (11)

d * [F.sup.2] = [J.sup.3] (12)

d[J.sup.3] = [0.sup.4] (13)

d[[epsilon].sup.1] = [F.sup.2] (14)

Equation (12) reveals that the electromagnetic charge J is a dual electromagnetic entity.

The projection of these continuous electromagnetic entities in Equations (11)-(14) on the cell complexes K and [??] provides the following discrete topological equations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (18)

where:

* [delta] is the coboundary operator.

* [[PHI].sup.2] = ([[PHI].sup.e], [[PHI].sup.b]) is the global electromagnetic flux corresponding to the electromagnetic strength 2-form [F.sup.2] = (E, B).

* [[??].sup.3] = ([[??].sup.[rho]], [[??].sup.j]) is the global electromagnetic charge corresponding to the electromagnetic charge 3-form J = ([rho], J).

* [U.sup.1] = ([U.sup.v], [U.sup.a]) is the global electromagnetic potential corresponding to the electromagnetic potential 1-form [[epsilon].sup.1] = (V, A).

* [[??].sup.2] = ([[??].sup.d], [[??].sup.h]) is the global electromagnetic charge-current potential corresponding to the dual of the electromagnetic strength 2-form *[F.sup.2] = (H, D).

We surmount by a tilde all global entities corresponding to dual electromagnetic entities. On their associated geometric objects endowed with inner/outer orientation, we reformulate the discrete equations as:

<[c.sup.i.sub.3], [delta][[PHI].sup.2]> = <[c.sup.i.sub.3], [0.sup.3]> (19)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

<[c.sup.i.sub.2], [delta][U.sup.1]> = <[c.sup.i.sub.2], [[PHI].sup.2]> (22)

where [c.sup.i.sub.p] ([[??].sup.i.sub.p]) designates a p-dimensional geometric object of index i endowed with inner (outer) orientation.

This formulation has important consequences [23] since the objective is to find a consistent discretization scheme for Maxwell's equations. Otherwise, these continuous and discrete equations permit to redraw the space-time classification diagrams, also called the Tonti diagrams [7-10, 24], as shown in the Figure 1.

By decoupling the electromagnetic fluxes to the electric and the magnetic fields, Equations (15) and (16) can be rewritten in a detailed manner:

[[PHI].sup.e]([partial derivative]S x I) + [[PHI].sup.b](S x [partial derivative]I) = [0.sup.3](S x I) (23)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

These equations are visualized graphically on the Figure 2. The Maxwell equations are given in a compact and elegant form [30-32] in terms of the four-dimensional differential form representation. This formalism appears to be an ideal tool for electromagnetic analysis [31], particularly for numerical algorithms. Moreover, Its geometric interpretation holds many pedagogical interests since it enables conjointly an intuitive and a dynamic visualization.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

2.4. Discretization of Metric Equations

The discrete topological equations were carried out without recourse to any approximation. Thus, discretization errors come only from the discretization of metric relations [9,33]. At first, we show how to perform the discretization of continuous operators * to be discrete operators [*]. On a complex K and its dual [??], we define the discrete operators [[*.sub.[epsilon]]], [[*.sub.[mu]]] and [[*.sub.[sigma]]] as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (25)

For a consistent discretization, the discrete Hodge operators [*] must complete some conditions. For example, on the usual Euclidean space, they should be symmetric and positive definite operators. When explicit schemes are preferred, it is necessary to verify that at least one of these matrixes needs to be diagonal. For an orthogonal grid mesh, we accomplish these conditions. The Voronoi-Delaunay mesh, by its construction, is a good example.

However, this discretization does not mean that [[*].sub.i] is constant. It simply means that the discrete operator is only depending on the cell where it is defined. In this way, the parameter--the permittivity [epsilon] as example--is an electromagnetic property of the corresponding object. Thus, a member function of the object--function UpdateEpsilon (...): Number as example--can adjust the parameter at any moment t and at any angular frequency w.

2.5. The General Finite Formulation

The discrete topological Equations (15)-(17) and (18) have the form:

[delta][a.sup.p] = [b.sup.p+1] (26)

where [a.sup.p] and [b.sup.p+1] are called cochains [25,28]. Cochains therefore constitute the discrete representation of the electromagnetic fields. This Equation (26) is known as the discrete strong form of the continuous Equations (3) and (4). By using the discrete generalized Stocke's theorem ([dagger]):

<[c.sup.i.sub.p+1], [delta][a.sup.p]> = <[partial derivative][c.sup.i.sub.p+1], [a.sup.p]> [for all][c.sup.i.sub.p+1] (27)

We get the discrete weak form:

<[c.sup.i.sub.p+1], [a.sup.p]> = <[c.sup.i.sub.p+1], [b.sup.p+1]> [for all][c.sup.i.sub.p+1] (28)

The boundary [partial derivative][c.sup.i.sub.p+1] is by definition the collection of its faces:

[partial derivative][c.sup.i.sub.p+1] = [[n.sub.p].summation over (j=1)] [c.sup.i.sub.p+1], [c.sup.j.sub.p]] [c.sup.j.sub.p] (29)

where:

(i) [n.sub.p] is the number of p-dimensional objects in the complex K.

(ii) [[c.sup.i.sub.p+1], [c.sup.j.sub.p]] is the incidence number resulting from the induced orientation of the cell [c.sup.i.sub.p+1] on the cell [c.sup.j.sub.p]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The linearity of the chain-cochain pairing gives:

[[n.sub.p].summation over (j=1)] [[c.sup.i.sub.p+1], [c.sup.j.sub.p]]<[c.sup.j.sub.p], [a.sup.p]> = ([c.sub.p+1], [b.sup.p+1]> [for all][c.sub.p+1] (30)

Thus, we rewrite Equations (19) and (20) as follows to obtain the finite formulation of the Element-Oriented Model:

[[n.sub.2].summation over (j=1)] [[alpha].sup.i,j.sub.3,2] <[c.sup.j.sub.2], [[PHI].sup.2]> = 0 (31)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

where [n.sub.2] is the number of the 2-dimensional cells in the complex K, [[alpha].sup.i,j.sub.3,2] = [[c.sup.i.sub.3], [c.sup.j.sub.2]] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

These Equations (31) and (32) govern the behaviors of the objects [c.sup.i.sub.3] and [[??].sup.i.sub.3] in the space-time structure previously described.

Decoupling the electromagnetic fluxes to the electric and the magnetic fields (see Figure 2) gives:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (34)

where:

* [c.sup.n + 1/2.sub.2] ([c.sup.n - 1/2.sub.2]) denotes the 2-dimensional geometric object (S x [partial derivative]I) localized at the instant t = n + 2 (t = n - 2).

* [[??].sup.n+1.sub.2]([[??].sup.n.sub.2]) denotes the 2-dimensional geometric object ([??] x [partial derivative][??]) localized at the instant [??] = n + 1 ([??] = n).

Notice that the electric flux [[PHI].sup.e] is computed on the cell [c.sup.j.sub.2] = [partial derivative]S x I during the time interval I = [n - 1/2 ,n + 1/2]. Similarly, [[??].sup.h] is computed during the time interval [??] = [n, n + 1] (see Figure 2). This statement constitutes the main disagreement with the most FDTD schemes where the flux [[PHI].sup.e]([[??].sup.h]) is calculated implicitly at the time instant t = n ([??] = n + 1/2). So, if we restrict these fluxes at an instant time, these equations become only approximation relations.

3. A FLEXIBLE FDTD SCHEME OF THE ELEMENT-ORIENTED MODEL

The computing restriction mentioned in the previous section also can be explained by considering the fluxes [[PHI].sup.e] and [[??].sup.h] constants on their associated surfaces [partial derivative]S x I and [partial derivative][??] x [??]. In this section, we examine a first idea to correct this disagreement by taking these fluxes variables on their surfaces. In other hand, the most known standard S-FDTD and nonstandard NS-FDTD methods use the Taylor form to approximate the differential operator acting on the electric and the magnetic fields. However, it is the solution that is sought, and hence it is redundant to approximate any function other than the solution itself [37, 38]. Unlike these methods, we attempt to incorporate any reasonable local approximating functions (such as plane wave, harmonic polynomial, sinusoidal, exponential, cylindrical or spherical harmonics, etc.) into the previous formulation of the Element-Oriented Model. In fact, solutions of many physical problems have salient local features that are qualitatively known a priori. In this situation, the accuracy of finite-difference methods in electromagnetic can be improved by approximating the solution itself and not the differential equations [37, 38].

For simplicity, we choose to develop our study (but not necessarily) on dual rectangular Cartesian grids. We focus our interest on the hypervolumetric primary cell H located between the coordinates (x = i, y = j, z = k, t = [t.sup.n - 1/2]) and (x = i + y = j + 1, z = k+1, t = [t.sup.n + 1/2]). It first base face is a volumetric primary cell ([V.sub.k] = [partial derivative]H[|.sub.z=k]) located at z = k between the coordinates (x = i, y = j, t = [t.sup.n - 1/2]) and (x = i + 1, y = j + 1, t = [t.sup.n + 1/2]) as illustrated on the Figure 3. The main basic idea is to consider a local solution [psi](x, y, z, t) defined on a cell. So, we assert that each component of the electric and magnetic fields can be written as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (35)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (36)

where [E.sub.[alpha]] and [B.sub.[alpha]] are the local magnitudes of the electric component [E.sub.[alpha]] and the magnetic induction component [B.sub.[alpha]] in a given direction [alpha] = x, y, z. These components [E.sub.[alpha]] and [B.sub.[alpha]] are defined at the center of the dual cells because they are deduced from the electric induction [??] and the magnetic field [??] also associated with the same cells. In general, at a node ([x.sub.i], [y.sub.j], [z.sub.k], [t.sup.n]), we have: [E.sub.[alpha]] = [f.sup.-1.sub.[epsilon]][??] and [B.sub.[alpha]] = [f.sub.[mu]][??]. We emphasize that the [psi] functions are defined on the dual cells which constitute also their supports.

[FIGURE 3 OMITTED]

From the Equations (33) and (34), we can deduce an algebraic set of equations. As example of demonstration, we choose ([double dagger]) to compute the flux [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] on the volumetric primary cell [V.sub.k]. On each face V of H and its dual, we obtain two algebraic equations [[PSI].sup.b] = f([[PSI].sup.e]) and [[??].sup.d] = f([[??].sup.h]).

On the left face of the electromagnetic object (see Figure 4(a)), The flux [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The flux is [[PSI].sub.e.sub.L!] is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (37)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

In the same manner we obtain:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (38)

Hence, we deduce the total flux [[PSI].sup.e.sub.L] on the left face:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (39)

On the right face, we obtain:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (40)

Therefore, the total flux [[PSI].sup.e.sub.R-L] = [[PHI].sup.e.sub.R] - [[PHI].sup.e.sub.L] is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (41)

The total flux [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (42)

The total flux of [B.sub.z] on the back-front faces is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (43)

Hence, we get an explicit form of the Equation (33) for the component [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

In the same manner, we obtain the explicit forms for the components [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] on the volumetric primary cells [V.sub.j] = [partial derivative]H[|.sub.y=j] and [V.sub.i] = [partial derivative]H[|.sub.x=i]. The components [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be easily computed on the dual cells [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] or simply deduced by permutation and translation.

Finally, we obtain a set of six equations that governs the behavior of the hypervolumetric space-time object H.

[FIGURE 4 OMITTED]

4. FROM THE FLEXIBLE FDTD SCHEME OF THE ELEMENT-ORIENTED MODEL TO THE STANDARD S-FDTD METHOD

In this section, we show that we can derive, under some assumptions, the standard S-FDTD method from the flexible FDTD scheme of the Element-Oriented Model. Let us consider the following assumptions:

(i) we assume that the dual mesh is orthogonal and regular. All intervals in the directions x, y, z, t, are equals to [DELTA]x, [DELTA]y, [DELTA]z and [DELTA]t respectively.

(ii) we assume that the normalized functions p are constants and equal to unity in their cells where they are defined as shown in Figure 5. Therefore, all functions [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] defined in Equations (37)-(41) become:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

All functions [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] defined in Equation (42) become: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]Ax At/4. All functions [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] defined in Equation (43) become:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

(iii) We assume that the field [E.sub.y] remains uniform over the intervals ]i - 1/2, i + 1/2], ]i + 1/2, i + 3/2]. Thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(iv) A linear approximation leads us to write:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (44)

Equation (39) becomes:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (45)

With the same assumptions on the right face of the volumetric cell [V.sub.k], the Equation (41) becomes:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (46)

(v) Assuming that Ex is uniform over the intervals ]j - 1/2, j + 2], ]j + 1/2,j + 3/2] and using a linear approximation, the Equation (42) becomes:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (47)

(vi) We consider [B.sub.z] uniform over the intervals ]n - 1 , n], ]n, n + 1]. Equation (43) becomes:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (48)

Equation (33) for the component [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] becomes:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (49)

Thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (50)

This formulation represents the discretization of the Faraday's law -[[partial derivative].sub.t]B = [nabla] x E by the conventional Yee-FDTD algorithm established on dual Yee-cells. The discrete derivative operator [d.sub.r] in a direction r is well:

[d.sub.r] x f = f(r + [DELTA]r) - f(r)/[DELTA]r

(vii) Assuming that the permittivity [epsilon] and the permeability [mu] are constants on every line c or [??] of any volumetric cell V or [??], we can write for any direction r = x, y, z at any point pt:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The operators [[*.sub.[epsilon]]] and [[*.sub.[mu]]] become:

[[*.sub.[epsilon]]] = [[[epsilon]].sub.3]; [[*.sub.[mu]]] = [[[mu]].sub.3]

where [[[epsilon]].sub.3] and [[[mu]].sub.3] are three order diagonal matrixes.

Equation (50) becomes:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (51)

Equation (34) (in the case of Q = 0) for the component [D.sub.z] becomes:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (52)

In the same manner, we find all remaining equations for the components [H.sub.x], [H.sub.y], [D.sub.x] and [D.sub.y] as previously mentioned.

Finally, the demonstration of the passage from the Element-Oriented Model to the standard S-FDTD is complete.

[FIGURE 5 OMITTED]

5. FROM THE FLEXIBLE FDTD SCHEME OF THE ELEMENT-ORIENTED MODEL TO THE NONSTANDARD NS-FDTD METHOD

To show the passage from the Element-Oriented Model to the nonstandard NS-FDTD, we maintain all previous assumptions in Section 4 except the assumption (ii). Therefore, we propose to consider the [psi] functions dependent on the directions of the cells where they are defined:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

These [psi] functions are dependent on the mesh density [N.sub.[lambda]] as shown in Figures 6 and 7. In the end, for a very fine mesh ([N.sub.[lambda]] [right arrow] +[infinity]). All [psi] functions will be equal to unity ([psi] [right arrow] 1). Hence, this formulation will be identical to the standard S-FDTD method. The flux of [E.sub.y] on the left face of the primary spatial temporal cell will be:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (53)

By computing and grouping fluxes:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (54)

we obtain:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (55)

It is clear that the discrete derivative operator has an unconventional form:

[d.sub.x]f(x) = f (x + h) - f(x)/[phi](h) where [phi](h) [right arrow] h when h [right arrow] 0 (56)

The same procedure can be applied to the other components to formulate completely this scheme. This formulation is well known as the nonstandard finite-difference time-domain method NS-FDTD [16, 21, 34]. The main rules to construct a nonstandard scheme are [16]:

* The order of derivative discretizations should be the same as the corresponding order derivatives in differential equations.

* The discrete form of the derivative must have a denominator in general a non-trivial function of the step discretization.

dy/dt [right arrow] [y.sub.k+1] - [psi](h)[y.sub.k]/[PHI](h)

where [PHI](h) = h + O([h.sup.2]) and [psi](h) = 1 + O(h)

* Nonlinear terms must be in general replaced by non-local representation: [x.sup.2] [right arrow] [x.sub.k+1][x.sub.k].

* The special conditions related to differential equations and/or their solutions are maintained for discrete formulations and/or their solutions, such as, for example, the time invariance t [right arrow] -t.

In this section, we have shown the passage from the Element-Oriented Model to the standard S-FDTD method.

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

6. A NEW NONSTANDARD FINITE-DIFFERENCE TIME-DOMAIN SCHEME

The main purpose from this section is to show the capability of the Element-Oriented Model to generate new schemes. In fact, we propose a new nonstandard NS-FDTD more accurate than the standard S-FDTD and the known nonstandard NS-FDTD methods. To formulate this new scheme, we will consider the same assumptions in Section 5. But, unlike the assumption (vi) in Section 4, we take the field [B.sub.z] variable on the interval ]n,n + 1]. We also take [D.sub.z] variable on the interval ]n + 1/2, n + 3/2].

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Near the instant t = n + 1/2, we take [B.sub.z][|.sup.n + 1/2] equal to its temporal mean value:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [tau] [member of] [0, +1].

In the same manner, near the instant t = n + 1, Dz[|.sup.n+1] becomes:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We can see that [B.sup.n + 1/2.sub.z] and [D.sup.n+1.sub.z] are corrected by the same function [[phi].sub.c]([DELTA]t, [tau]):

[[phi].sub.c]([DELTA]t, [tau]) = 2/[omega][tau][DELTA]t x sin([omega][tau][DELTA]t/2)

where:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

This function [[phi].sub.c]([DELTA]t, [tau]) is illustrated on the Figure 8 for some temporal steps.

[FIGURE 8 OMITTED]

Equation (52) becomes:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (57)

The new nonstandard temporal derivative operator is formulated as:

[d.sub.t] x f(t) = f(t + [DELTA]t) x [[phi].sub.c]([DELTA]t, [tau]) - f(t)/[phi]([DELTA]t) (58)

where: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This algorithm uses the new temporal derivative operator [d.sub.t] in Equation (58) and the spatial derivative operator [d.sub.x] in Equation (56).

Only these two operators [d.sub.t] and [d.sub.x] are sufficient to prove that the new nonstandard scheme generalizes both the standard and the nonstandard FDTD methods.

[FIGURE 9 OMITTED]

7. NUMERICAL EXAMPLES

A known severe test of any algorithm is Mie scattering [16]. In the Mie regime, the feature size of the scatterer is comparable to the wavelength. An infinite plane wave is incident from the left upon a perfectly conducting circular cylinder (Figure 9). In this case of cylindrical scatterers, the exact solution is given analytically by the Mie series. Using this analytic solution we are able to compute exact errors of any algorithm. The total electric field is calculated with the three methods (The standard S-FDTD, the nonstandard NS-FDTD and the new nonstandard NS-FDTD methods) and is compared with the analytic solution.

The exact analytical solution to the problem is given by [36]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (59)

where [E.sub.z] is the total electric field, [E.sub.0] is the amplitude of the plane wave, [k.sub.0] is the propagation constant in free space, [rho] is the radial distance from the center of the cylinder to the observation point, [theta] is the corresponding angle measured from the positive x-axis, [J.sub.n] is the Bessel function of the first kind, and [H.sup.(2).sub.n] is the Hankel function of the second kind.

[FIGURE 10 OMITTED]

[FIGURE 11 OMITTED]

Figure 10 shows the distribution of the total electric field Ez computed by the three numerical methods and compared with the exact solution.

Figure 11 shows the numerical errors based on the L1 and L2 norms. This proves that the new NS-FDTD deduced from the Element Oriented Model is more accurate than the standard and nonstandard FDTD methods with about the same computational cost as shown on the Figure 12.

[FIGURE 12 OMITTED]

[FIGURE 13 OMITTED]

In other hand, since our model is oriented object, we reuse the same program to develop several applications benefiting from the instantiation and the inheritance concepts. To illustrate the importance of this Element-Oriented Model, we present here another example of a photonic waveguide with three ports. This application is simply developed by instantiating some cylinders from the cylinder object of the previous problem. Figure 13 shows a simulation of this waveguide.

8. CONCLUSION

Based on a topological approach, we formulated a new Element-Oriented Model that conserves many physical informations. On dual staggered grids, we developed this model to obtain a flexible Finite-Difference Time-Domain using a general approximating function. We demonstrated that this scheme generalizes both the standard S-FDTD and the nonstandard NS-FDTD methods. Moreover, to confirm the productivity of this model, we presented a new nonstandard scheme. Its accuracy was proven by a numerical problem of a plane wave striking a perfectly conducting circular cylinder. Furthermore, the importance of this Element-Oriented Model was elucidated by instantiation of a cylinder object to produce a simulation of a photonic waveguide. We are persuaded that the Element-Oriented Model can also generate other new schemes. In future work, we will present a second new nonstandard NS-FDTD algorithm which is more accurate.

ACKNOWLEDGMENT

The authors would like to thank the anonymous reviewers for their helpful comments.

REFERENCES

[1.] Sadiku, M. N. O., Numerical Techniques in Electromagnetics, 2nd edition, CRC Press, 2001.

[2.] Salon, S. and M. V. K. Chari, Numerical Methods in Electromagnetism, Academic Press, 1999.

[3.] Rjasanow, S. and O. Steinbach, The Fast Solution of Boundary Integral Equations, Springer, 2007.

[4.] Baranger, J., J. F. Maitre, and F. Oudin, "Connection between finite volume and mixed finite element methods," RAIRO, Modelisation Math. Anal. Numer., Vol. 30, 445-465, 1996.

[5.] De La Bourdonnay, A. and S. Lala, "Duality between finite elements and finite volumes and Hodge operator," Numerical Methods in Engineering '96, 557-561, Wiley & Sons, Paris, 1996.

[6.] Bossavit, A. and L. Kettunen, "Yee-like schemes on staggered cellular grids: A synthesis between FIT and FEM approaches," IEEE Trans. Magn. , Vol. 36, No. 4, 861-867, 2000.

[7.] Teixeira, F. L., "Geometric aspects of the simplicial discretization of Maxwell's equations," Progress In Electromagnetics Research, Vol. 32, 171-188, 2001.

[8.] Tonti, E., "Finite formulation of the electromagnetic field," Progress In Electromagnetics Research, Vol. 32, 1-44, 2001.

[9.] Mattiussi, C., "An analysis of finite volume, finite element, and finite difference methods using some concepts from algebraic topology," J. Comp. Phys., Vol. 9, 295-319, 1997.

[10.] Gross, P. W. and P. R. Kotiuga, Electromagnetic Theory and Computation: A Topological Approach, Cambridge University Press, Cambridge, 2004.

[11.] Maxwell, J. C., A Treaties on Electricity and Magnetism, Clarendon Press, Oxford, 1892, Reprinted in 2002.

[12.] Stern, A., Y. Tong, M. Desbrun, and J. E. Marsden, "Computational electromagnetism with variational integrators and discrete differential forms," arXiv: 0707.4470 [math.NA], 2007.

[13.] Hiptmair, R., "Discrete Hodge-operators: An algebraic perspective," Progress In Electromagnetics Research, Vol. 32, 247-269, 2001.

[14.] Hiptmair, R., "Discrete Hodge operators," Numer. Math., Vol. 90, No. 2, 65-289, 2001.

[15.] Auchmann, B. and S. Kurz, "A geometrically defined discrete Hodge operator on simplicial cells," IEEE Trans. Magn., Vol. 42, No. 4, 643-646, 2006.

[16.] Mickens, R. E., Applications of Nonstandard Finite Difference Schemes, World Scientific, Singapore, 2000.

[17.] Marrone, M., "Computational aspects of the cell method in electrodynamics," Progress In Electromagnetics Research, Vol. 32, 317-356, 2001.

[18.] Garcia, S. G. and T.-W. Lee, "On the accuracy of the ADI-FDTD method," IEEE Antennas and Wireless Propagation Letters, Vol. 1, No. 1, 2002.

[19.] Ahmed, I., E. K. Chua, E. P. Li, and Z. Chen, "Development of the three dimensional unconditionally stable LOD-FDTD method," IEEE Trans. Antennas Propag., Vol. 56, No. 11, 3596-3600, 2008.

[20.] Zheng, F., Z. Chen, and J. Zhang, "A finite-difference time-domain method without the Courant stability conditions," IEEE Microw. Guided Wave Lett., Vol. 9, No. 11, 441-443, 1999.

[21.] Anguelov, R. and S. Lubuma, "On non-standard finite difference models of reaction-diffusion equations," Journal of Applied Mathematics, Vol. 175, No. 1, 2005.

[22.] Bossavit, A., "'Generalized finite differences' in computational electromagnetics," Progress In Electromagnetics Research, Vol. 32, 45-64, 2001.

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

[24.] Alotto, P., F. Freschi, and M. Repetto, "Multiphysics problems via the cell method: The role of Tonti diagrams," IEEE Trans. Magn., 2959-2962, Aug. 2010.

[25.] Tarhasaari, T., L. Kettunen, and A. Bossavit, "Some realizations of a discrete Hodge operator: A reinterpretation of finite element techniques," IEEE Trans. Magn., Vol. 35, No. 3, 1494-1497, 1999.

[26.] Tarhasaari, T. and L. Kettunen, "Topological approach to computational electromagnetism," Progress In Electromagnetics Research, Vol. 32, 189-206, 2001.

[27.] Truesdell, C. and R. A. Toupin, The Classical Field Theories, Harrdbuch der Physik, edited by S. Flugge, Vol. 311, 226-793, Springer-Verlag, Berlin, 1960.

[28.] Tarhasaari, T. and L. Kettunen, "Topological approach to computational electromagnetism," Progress In Electromagnetic Research, Vol. 32, 189-206, 2001.

[29.] Kirawanich, P., et al., "Methodology for interference analysis using electromagnetic topology techniques," Applied Physics Letters, Vol. 84, 2004.

[30.] Lindel, I. V., Differential Forms in Electromagnetics, IEEE Press, 2004.

[31.] Lindel, I. V., "Electromagnetic wave equation in differential-form representation," Progress In Electromagnetics Research, Vol. 54, 321-333, 2005.

[32.] Lindell, I. V., "Electromagnetic fields in self-dual media in differential-form representation," Progress In Electromagnetics Research, Vol. 58, 319-333, 2006.

[33.] Mattiussi, C., "The geometry of time-stepping," Progress In Electromagnetics Research, Vol. 32, 123-149, 2001.

[34.] Anguelov, R. and S. Lubuma, "Nonstandard finite-difference methods by nonlocal approwimations," Mathematics and Computer in Simulation, 2003.

[35.] Magrez, H. and A. Ziyyat, "Modelisation orientee objet en electromagntisme," Congres Mediterranen des Telecommunications CMT, Casablanca, 2010.

[36.] Balanis, C. A., Advanced Engineering Electromagnetics, Wiley, New York, 1989.

[37.] Pinheiro, H., J. P. Webb, and I. Tsukerman, "Flexible local approximation models for wave scattering in photonic crystal devices," IEEE Trans. Magn., Vol. 43, No. 4, 1321-1324, 2007.

[38.] Tsukerman, I., "A class of difference schemes with flexible local approximation," The Journal of Computational Physics, Vol. 211, No. 2, 659-699, 2006.

([dagger]) The continuous generalized Stocke's theorem is: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

([double dagger]) Since our mesh and our equations present many symmetries, only one equation is sufficient to deduce all remain equations by permutation and translation.

H. Magrez * and A. Ziyyat

Electronic and Systems Laboratory, Faculty of Sciences, Mohammed

1th University, Oujda 60000, Morocco

* Corresponding author: Hamid Magrez (magrezhamid@hotmail.com).

Received 1 May 2011, Accepted 3 November 2011, Scheduled 10 November 2011