# WHITNEY/NÉDÉLEC ELEMENTS METHOD APPROACH APPLIED IN THE MAXWELL'S EQUATIONS

Jean Eduardo Sebold jeansebold@ufpr.br/jean.sebold@ifc-araquari.edu.br <fnref n="1" />,<fnref n="2" />

Saulo Pomponet Oliveira saulopo@ufpr.br <fnref n="2" />

Luiz Alkimin de Lacerda alkimin@lactec.gov.br <fnref n="3" />

José Antonio Marques Carrer carrer@ufpr.br <fnref n="2" />

Rod. BR 280-km 27-CEP:89245-000-Araquari-SC-Brazil

Av. Cel. Francisco H. dos Santos, n/n - Jardim das Américas-CEP:81530-900-Curitiba-PR-Brazil

Av. Cel. Francisco H. dos Santos, n/n - Jardim das Américas-CEP:81530-900-Curitiba-PR-Brazil

**Abstract.** This work concerns Whitney and Nédélec finite element methods for time-harmonic Maxwell's equations. We review the derivation of the harmonic equations from full Maxwell's equations as well as their variational formulation, and build the Whitney and Nédélec element spaces, whose functions have continuous tangential components along the interface of adjacent elements. We study the dispersive behaviour of first-order Nédélec elements in two and three dimensions, in terms of the time frequency and the mesh element size, and present an explicit form for the discrete dispersion relation. Numerical experiments validate the performance of Whitney elements and Nédélec first order in a two-dimensional domain, that also illustrates the dispersion of the approximate solution with respect to the exact solution. The discrete dispersion relation for elements of the first order, show, through numerical evidence that the numerical phase velocity can be used as an error estimator in the Whitney and Nédélec finite element approximation, and thus, display an initial parameter h to the mesh refinement.

**Palavras-chave:** Whitney Elements, Nédélec Elements, Legendre hierarchical basis functions

## 1. INTRODUCTION

When in a numerical formulation is necessary and directly represent a generic form discretized vector quantity, the alternative feature that when using only basic functions of the nodal type is the separate treatment of each component of the field in question, which is individually simply reduce to scalar functions.

Arises from this fact that the first difficulties and concerns the continuity of the discretized vector quantity between adjacent elements of the finite element mesh. When you have a face in common, two neighboring and adjacent elements also share the nodes of the finite element mesh located on this face. Because these nodes belong simultaneously to two elements and because the approximation of vector quantities with nodal basis functions be performed component by component, it appears that the use of this approach implies the continuity of all components of the vector quantity in question in face shared by these neighboring elements. The final result of this is the fact that if perhaps each of the elements belonging to a material medium of different composition, this imposing continuity of all components produces, of course, a physically incorrect situation.

Such solutions are not guaranteed by, with the use of the nodal basis functions, the continuity of the derivatives of the interpolation functions between the elements of the mesh.

Fortunately, over time, way to overcome these difficulties have been devised. Among them, we highlight the development and improvement of techniques characterized by the use of another variety of finite element as an alternative to the original nodal approaches. Such alternatives have emerged in the work of Whitney (1957) and Nédélec (1980).

<footnote n="1">¹ Catarinense Federal Institute</footnote>

<footnote n="2">²Federal University of Paraná</footnote>

<footnote n="3">³Institute of Technology for Development</footnote>
Despite having used sets of edge vectors in a completely different context to the finite element method developed in this work, Whitney was one of the first to use a vector space of polynomials to generate such sets Monk (2003). These edge vectors zero order elements, in the context of the finite element method is an approximation order of elements that have constant tangential components on edges. In the literature we often find references to these elements as elements Whitney, Monk (2003).

A little later, Jean-Claude Nédélec presented some families of nonconforming finite elements in $\mathbb{R}^3$. We will see later that one of these families of finite elements is conforming elements in space $H(\text{curl}, \Omega) = \{u \in (L_2(\Omega))^2; \nabla \times u \in (L_2(\Omega))^2\}$. From these observations some applications of these elements in the approximation of Maxwell's equations and the elasticity equations were presented in Nédélec (1980). This innovative range of finite element is known as vector finite elements, or even elements of Nédélec.

Several authors have considered the dispersive behaviour of the finite element method, among which the most relevant to this work are: Christon (1999) considered that the dispersive behaviour of a variety of methods for finite wave equation of second order elements, and presented numerical comparisons between the discrete phase and group velocities of the exact values. Monk and Parrot (1994) considered the dispersive behaviour of the finite element first-order in triangular elements for Maxwell's equations when a refinement is made in the mesh. Monk and Cohen (1998) conducted a dispersion analysis of elements type Nédélec to time-dependent Maxwell's equations using a lumping together with the mass matrix of the tensor product of two and three-dimensional meshes. Ihlenburg and Babuška (1997) studied the dispersive properties of higher order finite element for the Helmholtz equation in one dimension, and obtained estimates for the approximation method to fifth order in which $\omega h < 1$. Numerical evidence has been presented, which led to the conjecture that the elements of order $p$ provide an approximation order $2p$ of dispersion relation when the mesh size $h$ tends to zero. Monk (2003) presented evidence, based on duality, for convergence of finite elements of Nédélec applied in a cavity problem for Maxwell's equations. This cavity was assumed to be a Lipschitz polyhedron, and the mesh was considered rectangular, but not uniform. Monk (1991) studied the use of Nédélec finite elements in $H(\text{curl}, \Omega)$ in the approximation of time harmonic Maxwell equations over a limited domain. This analysis revealed a rather complicated by the fact that the bilinear form was not considered coercive. This difficulty was circumvented by the use of discrete Helmholtz decomposition of the error vector. Ainsworth (2003) demonstrated that the numerical dispersion displays three different types of behavior, depending on the order of the method in relation to the mesh size and the wave number. These behaviours are described in the following sequence: Oscillation Phase, Transition Zone and Super-Exponential Decay. Ainsworth and Coyle (2001) studied a hierarchical basis functions set for the Galerkin descritization space $H(\text{curl}, \Omega)$ for both hybrid meshes containing quadrilaterals and triangles with a non-uniform arbitrary polynomial order. Ainsworth (2004) presented an argument showing that the discrete dispersion relation can be expressed in terms of an approximation of the scalar Helmholtz equation in one dimension. Moreover, clarified the relationship discrete dispersion in a valid dimension for arbitrary orders of approximation.

## 2. VECTORIAL FINITE ELEMENTS

The Sobolev spaces $H(\text{curl}, \Omega)$ plays a central role in the variational theory of Maxwell's equations, because according to Monk (2003) this space is the space of finite energy solutions, and this way we can guarantee the existence, uniqueness and regularity of physically significant discrete solutions Greenleaf, Kurylev, Lassas and Uhlmann (2007). Thus, it is appropriate to take this finite element space for a class of suitable finite elements subspace for the Maxwell's equations system. Another feature of this space is the choice of the finite element discretization, which is necessary for the tangential components of the field $E$ are continuous through the interface element, moreover, there is no obligation to the components of the normal are continuous.

Vectorial finite element can be used in complex geometries and also the presence of discontinuous electromagnetic properties. In the case of Maxwell's equations, the electric permittivity $\epsilon$ is discontinuous across the surface of a domain $\Omega$ of $\mathbb{R}^3$, it is known from electromagnetic theory that the tangential component of electric field $E$ is continuous across this surface, thus it is necessary that the tangential component of the field approximation $E_h$ is also continuous.

Choosing the Nédélec elements which guarantees an approximation $H(\text{curl})$-conforming, we see that the tangential component of the field approximation $E_h$ is continuous on the surface, in the case where two mesh elements having different material properties. In addition, the edge finite element have many interesting and challenging mathematical properties. Such properties are deeply explored in Monk (2003); Ainsworth (2003, 2004); Jin (2002); Monk (2003) among others.

We will restrict our study to the Whitney elements and Ndlec elements of the first order and applied in a convenient way, aiming to provide an initial understanding for readers interested in the foundations of this theory.

### 2.1. Degrees of freedom and Functions basis

We will give, from now on, a special attention to the two-dimensional finite elements. We define the space of polynomials $\hat{P}$ associated with the Nédélec element of order $p$ on the reference element $\hat{K} = (-1, 1)^2$ as

$$ \hat{P} = \mathbb{E}_p = \{ \hat{E} = (\hat{E}_1, \hat{E}_2); E_1 = \operatorname{span}\{\mathbb{Q}_{p,p+1}\} \quad \operatorname{e} \quad E_2 = \operatorname{span}\{\mathbb{Q}_{p+1,p}\} \}, \quad (1) $$
where $\mathbb{Q}_{p,q}$ is the monomials set of degree less than or equal to $p$ in $\hat{x}$ and of degree less than or equal to $q$ in $\hat{y}$, i.e.

$$ \mathbb{Q}_{p,q} = \{\hat{x}^i \hat{y}^j : 0 \le i \le p; 0 \le j \le q\} \qquad (2) $$

Let $\{L_k\}_{k=0}^p$ be the Legendre polynomials set, which are given by the Rodrigues formula

$$ L_k(\xi) = \frac{1}{2^k k!} \frac{d^{(k)}}{d\xi^{(k)}} [(\xi^2 - 1)^k] \quad 0 \le k \le p \qquad (3) $$

and, furthermore, we also define the set $\{l_k\}_{k=0}^{p+1}$ as

$$ \left\{ \begin{array}{ll} l_k(\xi) = \frac{1}{2}(1 + \xi_k \xi), & k = 0, 1, \\ l_k(\xi) = \int_{-1}^{\xi} L_{k-2}(t) dt, & k = 2, \dots, p+1 \end{array} \right. \qquad (4) $$

where $\xi_0 = -1$, $\xi_1 = 1$ and $l_k(-1) = l_k(1) = 0$ for $k > 1$. In fact, for $\xi = -1$ the statement is clear. For $\xi = 1$, $l_k(1) = 0$ by the orthogonality of Legendre polynomials. Set $\{l_k\}_{k=0}^{p+1}$ is used for setting hierarchical basis function, which are useful for adaptive finite element methods, Adjerid (2001).

We begin by defining the degrees of freedom on a reference element $\hat{K} \subset \mathbb{R}^2$, according to Nédélec (1980). In this work we will adopt quadrilateral elements, thus we consider $\hat{K} = (-1, 1)^2$.

One way to ensure continuity between two elements with the same interface, is to select the degrees of freedom on an edge $\gamma$ to be the weighted moments of the tangential component of a field $\hat{E} \in \mathbb{E}_p$. Degrees of freedom at the edges of the reference element is defined by linear functional $\alpha_{\gamma} \in \mathcal{P}'$

$$ \alpha_{\gamma}(\mathbf{u}) = \int_{\gamma} (\mathbf{t} \cdot \mathbf{u}) \phi \, ds \quad \forall \phi \in \mathbb{P}_{k-1}(\gamma) \quad k = 1, \dots, p+1 \qquad (5) $$

where $\mathbf{t}$ is the unit tangent vector to the edge, $\mathcal{P}'$ denotes the dual space $\mathcal{P}$ and $\mathbb{P}_{k-1}$ is the linear space of polynomials, the degree of which is less or equal to $k-1$. We have a total of $4(p+1)$ degrees of freedom in the edges. Degrees of freedom within the reference element are defined by the linear functional $\alpha_{int} \in \mathcal{P}'$

$$ \alpha_{int}(\mathbf{u}) = \int_{\hat{K}} (\mathbf{u} \cdot \phi) \, d\hat{K} \quad \forall \phi = \begin{bmatrix} \phi_1 \\ \phi_2 \end{bmatrix} \quad \text{with } \phi_1 \in \mathbb{Q}_{k-2,k-1} \text{ and } \phi_2 \in \mathbb{Q}_{k-1,k-2} \quad k = 2, \dots, p+1 \qquad (6) $$

We have a total of $2p(p+1)$ degrees of freedom within the reference element $\hat{K}$.

Using the functions $L_k$ and $l_k$ defined above, we define the basis functions associated to the degrees of freedom at the edges as

$$ \left\{ \begin{array}{ll} \phi_{i,j,1}^{\gamma}(\hat{x}) = L_i(\hat{x})l_j(\hat{y})e_1 & i = 0, \dots, p; \quad j = 0, 1 \\ \phi_{i,j,2}^{\gamma}(\hat{x}) = l_j(\hat{x})L_i(\hat{y})e_2 & \end{array} \right. \qquad (7) $$

where $e_1 = [1 \ 0]^T$ and $e_2 = [0 \ 1]^T$ denote the canonical basis of $\mathbb{R}^2$ and $\hat{x}, \hat{y}$ are the coordinates of $\hat{x}$. Basis functions associated to the degrees of freedom inside of element are complement functions on (7) when using $l_k$ with $k \ge 2$:

$$ \left\{ \begin{array}{ll} \phi_{i,j,1}^{int}(\hat{x}) = L_i(\hat{x})l_j(\hat{y})e_1 & i = 0, \dots, p; \quad j = 2, \dots, p+1, \\ \phi_{i,j,2}^{int}(\hat{x}) = l_j(\hat{x})L_i(\hat{y})e_2 & \end{array} \right. \qquad (8) $$

Despite the representation (7)-(8) have the convenience of be compact and general, we use a different notation for the functions of the edges; This rating takes into account the edges of the reference element chosen:

$$ \phi_i^{\gamma_1}(\hat{x}) = L_i(\hat{x})l_0(\hat{y})e_1; \quad \phi_i^{\gamma_2}(\hat{x}) = L_i(\hat{x})l_1(\hat{y})e_1; \qquad (9) $$

$$ \phi_i^{\gamma_3}(\hat{x}) = L_i(\hat{y})l_0(\hat{x})e_2; \quad \phi_i^{\gamma_4}(\hat{x}) = L_i(\hat{y})l_1(\hat{x})e_2; $$

Thus, the basis function associated with the edges satisfy

$$ \int_{\gamma_u} (\mathbf{t} \cdot \phi_i^{\gamma_j}) L_k ds = \delta_{lj} \delta_{ik} ||L_i||^2 \quad \text{with } j = 1, 2, 3, 4 \qquad (10) $$

Similarly, the basis functions associated to inside the element satisfies

$$ \int_{\hat{K}} L_m(\hat{x})l_n(\hat{y})\phi_{i,j,1}^{int} \cdot e_1 d\hat{x}d\hat{y} = \delta_{mi}\delta_{nj}||L_i||^2||l_j||^2 \qquad (11) \\ \int_{\hat{K}} L_m(\hat{y})l_n(\hat{x})\phi_{i,j,2}^{int} \cdot e_2 d\hat{x}d\hat{y} = \delta_{mi}\delta_{nj}||L_i||^2||l_j||^2 $$

A Nédélec element $(\mathcal{P}, K, \Sigma)$ on the physical domain $K$ is built from the reference element $(\hat{\mathcal{P}}, \hat{K}, \hat{\Sigma})$ as follows:
Let $E$ be a vector field. The $m^{th}$ degree of freedom on the edge $\zeta$ is given by the linear mapping
$$ \boldsymbol{E} \rightarrow \int_{\zeta} v_m \boldsymbol{E} \cdot d\boldsymbol{x}, $$

where the weight function $v_m$ is chosen to coincide with the $m^{th}$ Legendre polynomial when the edge $\zeta$ is parametrized by $t \in (-1, 1)$. In particular, let $\zeta$ be an oriented edge of an element in the mesh with endpoints $\boldsymbol{x}_{v_a}$ and $\boldsymbol{x}_{v_b}$, where the indices $v_a$ and $v_b$ correspond to adjacent vertices in the quadrilateral element. Now, we introduce the parameterization on the edge as

$$ \boldsymbol{x}(t) = (x(t), y(t)) = \frac{1}{2}(1-t)\boldsymbol{x}_{v_a} + \frac{1}{2}(1+t)\boldsymbol{x}_{v_b} \quad \text{with} \quad t \in (-1, 1) \qquad (12) $$

Hence, using equation 5.34 found in Kaplan (1970), we obtain

$$ \int_{\zeta} v_m \boldsymbol{E} \cdot d\boldsymbol{x} = \int_{-1}^{1} L_m(t) \left( E_1 \frac{dx}{dt} + E_2 \frac{dy}{dt} \right) dt = \int_{-1}^{1} L_m(t) \boldsymbol{E} \cdot \boldsymbol{\sigma}_{\zeta} dt, \qquad (13) $$

where $\boldsymbol{\sigma}_{\zeta} = (\frac{dx}{dt}, \frac{dy}{dt})$ is the tangent vector on the edge $\zeta$.

Let $F_K$ be a mapping of the reference element $K$, $F_K : \hat{K} \rightarrow K$, defined by

$$ F_K(\hat{\boldsymbol{x}}) = J_K \hat{\boldsymbol{x}} + \boldsymbol{b}, \qquad (14) $$

where $J_K$ is an invertible square matrix and $\boldsymbol{b}$ is a translation vector. Note that $F_K$ is a differentiable bijection, and the Jacobian matrix transformation is given by $dF_K = J_K$. It has from the covariant transformation, Kaplan (1970), that the tangent vectors $\boldsymbol{\sigma}$ and $\boldsymbol{\tau}$ are related by

$$ \boldsymbol{\sigma}(\boldsymbol{x}) = J_K(\hat{\boldsymbol{x}})\boldsymbol{\tau}(\boldsymbol{\xi}) $$

therefore, we can represent (13) as

$$ \int_{-1}^{1} L_m(t) \boldsymbol{E} \cdot (J_K \boldsymbol{\tau}) dt $$

and using the mapping $F_\gamma$, restricted to reference edge $\gamma$ and with $m = k$, we obtain

$$ \int_{-1}^{1} L_k(t) \boldsymbol{E} \cdot (J_K \boldsymbol{\tau}) dt = \int_{-1}^{1} L_k(s) J_K^T \boldsymbol{E} \cdot \boldsymbol{\tau} ds $$

Note that the electric field $\hat{\boldsymbol{E}}$ on a reference element $\hat{K}$ is related to the electric field $\boldsymbol{E}$, defined on physical element $K$ by covariant transformation Kaplan (1970)

$$ \boldsymbol{E}(\boldsymbol{x})|_{K} = J_{K}^{-T} \hat{\boldsymbol{E}}(\hat{\boldsymbol{x}}), \qquad \boldsymbol{x} = F_{K}(\hat{\boldsymbol{x}}). \qquad (15) $$

Similarly, the global basis function $\phi_{i,j}^{(d)}$ ($d = 1, 2$) corresponding to the local basis function $\hat{\phi}_{i,j}^{(d)}$ on the reference element is defined by

$$ \phi_{i,j}^{(d)}(\boldsymbol{x})|_{K} = J_{K}^{-T} \hat{\phi}_{i,j}^{(d)}(\hat{\boldsymbol{x}}), \qquad \boldsymbol{x} = F_{K}(\hat{\boldsymbol{x}}). \qquad (16) $$

The degrees of freedom on the edges ensures that $(K, \mathcal{P}, \Sigma)$ is $H(\text{curl}, \Omega)$-conforming; complementing them with the degrees of freedom inside of element. We also ensure that $(K, \mathcal{P}, \Sigma)$ is unisolvente (see Theorem 5 of Nédélec (1986)).

## 2.2. Basis Function for Whitney Elements ($p = 0$)

An emphasis on the definition of $E_p$ for the case $p = 0$, implies $Q_{0,1} = \{1, \hat{y}\}$ and $Q_{1,0} = \{1, \hat{x}\}$, therefore

$$ E_0 = \{ \hat{\boldsymbol{E}} = (\hat{E}_1, \hat{E}_2); \hat{E}_1 = \operatorname{span}\{1, \hat{y}\} \text{ and } \hat{E}_2 = \operatorname{span}\{1, \hat{x}\} \} \qquad (17) $$

For Whitney elements there is only one degree of freedom at each edge, equation (5), thus we can define the basis functions as

$$ \begin{aligned} \phi_0^{\gamma_1}(\hat{x}) &= L_0(\hat{x})l_0(\hat{y})e_1; \quad \phi_0^{\gamma_2}(\hat{x}) = L_0(\hat{x})l_1(\hat{y})e_1; \\ \phi_0^{\gamma_3}(\hat{x}) &= L_0(\hat{y})l_0(\hat{x})e_2; \quad \phi_0^{\gamma_4}(\hat{x}) = L_0(\hat{y})l_1(\hat{x})e_2, \end{aligned} \qquad (18) $$

where $e_1$ and $e_2$ denote the canonical basis of $\mathbb{R}^2$. Zero-order functions can be characterized by the fact that they have free divergent, i.e., $\nabla \cdot \phi_0^{\gamma_i} = 0$, and also for having the constant tangential component on each edge $\gamma_i$, Jin (2002). Figure 1(a) shows in plane $\hat{x}\hat{y}$, the degrees of freedom distributed on the edges of the reference element.
## 2.3. Basis Function for Nédélec ($p = 1$)

Space of polynomials associated with the element of Nédélec of order $p = 1$ on the reference element $\hat{K} = (-1, 1)^2$ is defined as

$$ \mathbb{E}_1 = \{ \hat{\mathbf{E}} = (\hat{E}_1, \hat{E}_2); \hat{E}_1 = \operatorname{span}\{\mathbb{Q}_{1,2}\} \text{ and } \hat{E}_2 = \operatorname{span}\{\mathbb{Q}_{2,1}\} \}, \quad (19) $$

where

$$ \mathbb{Q}_{1,2} = \{1, \hat{y}, \hat{y}^2, \hat{x}, \hat{x}\hat{y}, \hat{x}\hat{y}^2\} \quad \text{e} \quad \mathbb{Q}_{2,1} = \{1, \hat{x}, \hat{x}^2, \hat{y}, \hat{y}\hat{x}, \hat{y}\hat{x}^2\} \quad (20) $$

$\mathbb{E}_1$ is constructed by increasing the space $\mathbb{E}_0$ – Hierarchical basis. Furthermore, we note that there are four new basis functions on edges that are generated by the elements $\hat{x}, \hat{x}\hat{y} \in \mathbb{Q}_{1,2}$ and $\hat{y}, \hat{y}\hat{x} \in \mathbb{Q}_{2,1}$, which will be added to the edges together with the basis functions for elements of Whitney.

Unlike the basis functions for Whitney elements, here, we have interiors basis functions generated by elements $\hat{x}, \hat{y}^2, \hat{x}\hat{y}^2 \in \mathbb{Q}_{1,2}$ and $\hat{y}, \hat{x}^2, \hat{y}\hat{x}^2 \in \mathbb{Q}_{2,1}$. Figure 1(b) shows the distribution of degrees of freedom on the quadrilateral element.

<figure><img src="image_1.png" /><figcaption>Figure 1. (a) Degrees of freedom for Whitney elements. The arrows indicate the zero-order moment of the tangential field $\hat{\mathbf{E}}$. Note that the direction of the tangent remains in the anti-clockwise around the element; (b) The arrows indicate the degrees of freedom of the elements of order $p = 1$ distributed over edges. The black squares in the interior of the element represent the degrees of freedom necessary for a quadrilateral.</figcaption></figure>

As we can see, the elements of order $p = 1$ are armed with degrees of freedom on the edges and inside the element. Thus, the basic functions of these elements are organized as follows:

1. Basis function in the edges as

$$ \left\{ \begin{array}{ll} \phi_0^{\gamma_1}(\hat{x}) = L_0(\hat{x})l_0(\hat{y})e_1; & \phi_0^{\gamma_2}(\hat{x}) = L_0(\hat{x})l_1(\hat{y})e_1; \\ \phi_0^{\gamma_3}(\hat{x}) = L_0(\hat{y})l_0(\hat{x})e_2; & \phi_0^{\gamma_4}(\hat{x}) = L_0(\hat{y})l_1(\hat{x})e_2 \\ \phi_1^{\gamma_1}(\hat{x}) = L_1(\hat{x})l_0(\hat{y})e_1; & \phi_1^{\gamma_2}(\hat{x}) = L_1(\hat{x})l_1(\hat{y})e_1; \\ \phi_1^{\gamma_3}(\hat{x}) = L_1(\hat{y})l_0(\hat{x})e_2; & \phi_1^{\gamma_4}(\hat{x}) = L_1(\hat{y})l_1(\hat{x})e_2 \end{array} \right. \quad (21) $$

Note that the first four basis functions of edge are the same basic functions for elements of order $p = 0$. The basis functions in (21) are built using the weight $v_k = L_k$, so satisfying (10).

2. Basis functions inside $\phi_{i,j,1}^{int}$ and $\phi_{i,j,2}^{int}$, with $i = 0, 1$ and $j = 2$, as

$$ \left\{ \begin{array}{ll} \phi_{0,2,1}^{int_x}(\hat{x}) = L_0(\hat{x})l_2(\hat{y})e_1; & \phi_{1,2,1}^{int_x}(\hat{x}) = L_1(\hat{x})l_2(\hat{y})e_1; \\ \phi_{0,2,2}^{int_y}(\hat{x}) = L_0(\hat{y})l_2(\hat{x})e_2; & \phi_{1,2,2}^{int_y}(\hat{x}) = L_1(\hat{y})l_2(\hat{x})e_2 \end{array} \right. \quad (22) $$

Basis functions in (22) are constructed to satisfy (11).
### 3. PERFORMANCE OF WHITNEY/NÉDÉLEC ELEMENTS

Now, in order to illustrate the performance of Whitney/Nédélec elements, consider a numerical example of a propagation of a plane wave through a square domain $\Omega = (0, 1)^2$ with boundary $\Gamma$. Let the fixed frequency $\omega$ and sufficiently smooth function $g$ be. Based on time harmonic Maxwell's second order system, Monk (2003), we approximate the electric field $\mathbf{E} \in H(\text{curl}, \Omega)$ such that

$$ (\nabla \times \mathbf{E}, \nabla \times \mathbf{v}) - \omega^2(\mathbf{E}, \mathbf{v}) = \mathbf{0} \qquad (23) $$

for all $\mathbf{v} \in H_0(\text{curl}, \Omega) = \{\mathbf{u} \in H(\text{curl}, \Omega); t \cdot \mathbf{u} = 0\}$, where $(\cdot, \cdot)$ defines the inner product in $(L_2(\Omega))^2$, Kreyszig (1978). Moreover, the operation $\nabla \times \mathbf{E}$ is the surface scalar rotational, Boffi and Perugia (1999); Monk (2003), defined for any vector function $\mathbf{u} = (u_1(x, y), u_2(x, y))$ as $\nabla \times \mathbf{u} = \frac{\partial u_2}{\partial x} - \frac{\partial u_1}{\partial y}$. Problem (23) is subject to the boundary condition

$$ t \cdot \mathbf{E} = g \quad \text{on } \Gamma, \qquad (24) $$

where $t$ is the unit tangent vector on the boundary $\Gamma$.

Let us choose the function $g$ so that we have the exact solution known. Following Ainsworth (2004), let $\tilde{\mathbf{E}} = i\nabla \times e^{i\xi \cdot \mathbf{x}}$ be a plane wave propagating in the direction $\xi = (4\pi, 4\pi)$. If we choose $\omega^2 = 32\pi^2$ and $g = t \cdot \tilde{\mathbf{E}}$, then $\tilde{\mathbf{E}}$ is the exact solution of the problem (23)-(24).

Consider that the domain $\Omega = (0, 1)^2$ is discretized on a uniform mesh of square elements of side $h$, Figure 2, and uniform $p$-order Nédélec elements are used to define the space of finite elements $V_{hp}$.

The finite element approximation is seeking $E_{hp} \in V_{hp}$ such that

$$ (\nabla \times E_{hp}, \nabla \times \mathbf{v}) - \omega^2(E_{hp}, \mathbf{v}) = \mathbf{0} \qquad (25) $$

for all $\mathbf{v} \in V_{hp} \cap H_0(\text{curl}, \Omega)$. Essential boundary conditions are imposed requiring that over all edges $\gamma \subset \Gamma$, we have

$$ \int_{\gamma} (t \cdot E_{hp} - g) \mathbf{v} \, ds = 0 \qquad (26) $$

for all $\mathbf{v} \in \mathbb{P}_p(\gamma)$, where $\mathbb{P}_p$ denotes the one-dimensional space of polynomials of degree at most $p$ in the arc length.

<figure><img src="image_2.png" /><figcaption>Figure 2. Example particular mesh used in the numerical experiments, with $n^2 = 16$, being $n$ the number of elements in the base domain, i.e. $h = 1/n$.</figcaption></figure>

Following, we will opt for making approximations to $\tilde{\mathbf{E}}$ using Whitney elements and first order Nédélec elements. Present some numerical experiments concerning the problem seen earlier in this section. Figures 3-5 highlights the approximation of the real part of the second component of the exact solution $\tilde{\mathbf{E}}$ to the problem (23)-(24).

<figure><img src="image_3.png" /><figcaption>Figure 3. (a) Exact solution; (b) Numerical solution with Whitney elements and $h = 1/16$; (c) Variation along the main diagonal of the Figures 3(a) and 3(b) of exact solution (dashed line) and numerical solution (continuous line), respectively.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 4. (a) Exact solution; (b) Numerical solution with Whitney elements and $h = 1/32$; (c) Variation along the main diagonal of the Figures 4(a) and 4(b) of exact solution (dashed line) and numerical solution (continuous line), respectively.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5. (a) Exact solution; (b) Numerical solution with first order Nédélec elements and $h = 1/4$; (c) Variation along the main diagonal of the Figures 5(a) and 5(b) of exact solution (dashed line) and of numerical solution (continuous line), respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6. (a) Exact solution; (b) Numerical solution with first order Nédélec elements and $h = 1/16$; (c) Variation along the main diagonal of the Figures 6(a) and 6(b) of exact solution (dashed line) and of numerical solution (continuous line), respectively.</figcaption></figure>

## 4. SELECTING THE MESH PARAMETER

In section earlier we saw that the error of approximation by finite elements may cause a phase difference with respect to the exact solution, see Figure 5(c). This effect depends not only on the *h* mesh parameter, but also the temporal frequency ω as illustrated in Figures 7(a) and 7(b). This fact can be analyzed by studying the dispersive properties of numerical solution in an infinity mesh.

Expressions that represent: (i) the dispersion relation for the equation of the electric field $\nabla \times \nabla \times E - \omega^2 E = 0$; (ii) the discrete dispersion relation elements of Nédélec of dimension *d* = 2 of arbitrary order (Whitney elements when
<figure><img src="image_7.png" /><figcaption>Figure 7. (a) Variation along the main diagonal ($h = 1/8$) of exact solution (dashed line) and of numerical solution with Nédélec elements of the first order (solid line) with $\omega^2 = 32\pi^2$ and with (b) $\omega^2 = 128\pi^2$.</figcaption></figure>

$p = 0$), are found in Ainsworth (2003, 2004), and are given respectively by

$$ \omega = |\xi| \qquad (27) $$

and

$$ \omega^2 = \sum_{i=1}^{d=2} \omega_{hp}(\xi_i)^2 \quad \text{with} \quad h \to 0 \qquad (28) $$

where $\xi_{1,2}$ are the components of the wave vector $\xi$.

The dispersion relation (27), tells us that the phase velocity of the continuous problem is $c = 1$. In the case of an approximation by Whitney elements, according to Ainsworth (2004), we can express each plot $\omega_{h0}(\xi_i)$, with $i = 1, 2$, in (28) as

$$ \omega_{h0}(\xi)^2 = \frac{6}{h^2} \left( \frac{1 - \cos(h\xi)}{2 + \cos(h\xi)} \right), \qquad (29) $$

while that for first order Nédélec elements one has

$$ \omega_{h1}(\xi)^2 = \frac{1}{h^2} \left( \frac{16\cos(\xi h) + 104 + \sqrt{(16\cos(\xi h) + 104)^2 - 4(\cos(\xi h) - 3)(\cos(\xi h) - 1)}}{6 - 2(\cos(\xi h))} \right) \qquad (30) $$

Setting the numerical phase velocity as

$$ C = \frac{\omega_{hp}(\xi)}{\xi} \qquad (31) $$

we can see in Figure 8(a) the comparison between the exact phase velocity of the continuous problem and the numerical phase velocity when we substitute the expressions (29) and (30), one and then another, in equation (31). Figure 8(b) is a closer view of Figure 8(a) that estimates, for example, the largest possible value of the parameter $h$ so that the estimated phase error is less than 0.01 %. To do so, simply observe the points where the velocity curves reach the value 1.0001, i.e., $\xi h \approx 0, 05$ and $\xi h \approx 0, 62$ for Whitney and Nédélec elements, respectively.

We will validate the estimates obtained for $h$ from the numerical experiments made in section 3. Use in this section $\xi = (4\pi, 4\pi)$, then $\xi = 4\pi$ and thus recommended for $h$ values are: $h \approx \frac{0.05}{4\pi} \approx \frac{1}{251}$ (Whitney elements) and $h \approx \frac{0.62}{4\pi} \approx \frac{1}{20}$ (Nédélec elements).

In fact, for example, from the parameter $h \le \frac{1}{20}$ phase difference in the experiment with Nédélec elements becomes negligible within the error suggested. This shows that the dispersive effects caused by numerical approximation by Whitney($p = 0$) and Nédélec($p = 1$) elements can be controlled according to the choice of initial parameter $h$ estimated by this analysis.

Let $\epsilon$ be the relative error in the phase velocity given by $|1 - C|$. Equations (29) and (30) are employed to show the phase velocity approximation for a fixed $\xi$ and an increasing $n$. This is shown in Figures 9(a) and 9(b) for three different $\xi$ values. As expected, it is clear that improving the approximation for $c = 1$ requires smaller $h$ values for larger frequency numbers, similar way the approximation by vector finite element, see Figures 7(a) and 7(b). Therefore, we can characterize the error in the phase velocity as error estimator in the vector finite element approximation.
<figure><img src="image_8.png" /><figcaption>Figure 8. (a) Numerical phase velocity for Whitney and Nédélec elements, and phase velocity $c = 1$.; (b) Closer view Figure 8(a)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9. (a) Error between the phase velocity of the continuous problem and the numerical phase velocity for Whitney elements; (b) Nédélec elements of first order. The number of elements in the mesh is given by $n^2$, where $h = 1/n$. (see Figure 2).</figcaption></figure>

## 5. ACKNOWLEDGMENTS

Acknowledgments in portuguese language: Agradecemos à Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - CAPES, ao Programa de Pós-Graduação em Métodos Numéricos em Engenharia - PPGMNE - UFPR, ao Instituto de Tecnologia para o Desenvolvimento - LACTEC, Instituto de Matemática Industrial - IMI - Curitiba/Paraná e ao Instituto Federal Catarinense - IFC pelo apoio.

## REFERENCES

* Adams, R. 1975. *Sobolev Spaces*. New York, USA: Academic Press. 300p.

* Ainsworth, M. 2003. Discrete dispersion for hp-version finite element approximation at high wavenumber. *SIAM J. Numer. Analysis*, vol.42, 553–575.

* Ainsworth, M. 2004. Dispersive properties of high - order Nédélec/edge element approximation of the time - harmonic Maxwell equations. *Philosophical transactions of the The Royal Society of London*, vol.362, 471–491.

* Ainsworth, M. and Coyle, J. 2001. Hierarchical hp-edge element families for Maxwell's equations on hybrid quadrilateral and triangular meshes. *Comput. Methods Appl. Mech. Engrg*, vol. 190, 6709–6733.

* Adjerid, S. 2001. Hierarchical finite element bases for triangular and tetrahedral elements. *Comput. Methods Appl. Mech. Engrg*, vol.190, 2925–2941.

* Boffi, D. and Perugia, I. 1999. Computation model of electromagnetic resonator: analysis of edge element approximation. *SIAM J. Numer. Ana*, vol. 36, 1264–1290.

* Brenner, S. and Scott, L. 1994. *The Mathematical Theory of Finite Element Methods*. New York, USA: Springer-Verlag. 380p.

* Brezzi, F. and Fortin, M. 1991. *Mixed and Hybrid Finite Element Methods*. New York, USA: Springer-Verlag. 421p.

* Christon, M. 1999. The influence of the mass Matrix on the dispersive nature of the semi-discrete, second-order wave equation. *Comput. Methods Appl. Mech. Engrg*, vol.173, 146–166.

* Ciarlet, G. 1978. *The Element Method for Elliptic Problems*. New York, USA: North-Holland Pub. 551p.
Colton, D. and Kress, R. 1983. *Integral Equations Methods in Scattering Theory*. New York, USA: John Wiley and Sons InC. 378p.

Girault, V. and Raviart, P. 1986. *Finite Element Methods for Navier-Stokes Equations - Theory and Algorithms*. Berlin, Deutschland: Springer-Verlag. 423p.

Girault, V., Bernardi, C., Amrouche, C. and Dauge, M. 1998. Vector potentials in three-dimensional non-smooth domains. *Mathematical Methods in the Applied Sciences*, vol. 21, 823-864.

Greenleaf, A. Kurylev, Y. Lassas, M. and Uhlmann, G. 2007. Full-Wave Invisibility of Active Devices at All Frequencies. *Communications in Mathematical Physics*, vol. 275, 749-789.

Grisvard, P. 1985 *Elliptic Problems in Nonsmooth Domains*. Massachusetts, USA: Pitman Publishing INC. 287p.

Haueisen, J. 1997. Influence of Tissue Resistivities on Neuromagnetic Fields and Electric Potentials Studied with a Finite Element Model of the Head *IEEE Transactions on Biomedical Engineering*, vol. 44, 727-734.

Ihlenburg, F. and Babuška, I. 1997. Finite element solution of the Helmholtz equation with high wave number. *Society for Industrial and Applied Mathematics*, vol.34, 315-358.

Jin, J. 2002. *The Finite Element Method in Electromagnetism Second Edition*. New York, USA: John Wiley and Sons InC. 389p.

Kaplan, W. 1970. *Advanced Calculus*. New York, USA: House of Electronics Industry. 754p.

Kreyszig, E. 1978. *Introductory Functional Analysis With Applications*. New York, USA: John Wiley and Sons. 704p.

Lima, E. 2006. *Um Curso de Análise Vol. 2*. Rio de Janeiro, Brazil: IMPA-Projeto Euclides. 546p.

Li, Y and Dai, S. 2011. Finite element modelling of marine controlled-source electromagnetic responses in two-dimensional dipping anisotropic conductivity structures. *Geophysical Journal International*, vol.185, 622–636.

Monk, P. 2003. *Finite Element Methods for Maxwell's Equations*. New York, USA: Oxford Science Publications. 465p.

Monk, P. 1991. A finite element method for approximating the time-harmonic Maxwell equation. *Numerische Mathematik*, vol.63, 243–261.

Monk, P. 2003. A simple Proof of Convergence for an Edge Element Discretization of Maxwell Equations. *Lecture Notes in Computational Science and Engineering*, vol.28, 127–141.

Monk, P. and Parrot, A. 1994. Dispersion analysis of finite element methods for Maxwell equations. *SIAM J. Sci. Comput.*, vol.15, 916–937.

Monk, P. and Cohen, G. 1998. Gauss point mass lumping schemes for Maxwell equations. *Numerical Methods for PDEs.*, vol.14, 63–88.

Nečas, J. 1983 *The Introduction to the Theory of Nonlinear Elliptic Equations*. New York, USA: John Wiley and Sons. 264p.

Nédélec, J. 1980 Mixed finite elements in $\mathbb{R}^3$. *Numerische Mathematik*, vol.35, 315-341.

Nédélec, J. 1986. A new family of mixed finite elements in $\mathbb{R}^3$. *Numerische Mathematik*, vol.50, 57–81.

Reddy, B. 1986. *Functional Analysis and Boundary-Value Problems: an Introductory Treatment*. New York, USA: John Wiley and Sons. 426p.

Smith, J. 2003. Human exposure assessment in the near field of GSM base-station antennas using a hybrid finite element/method of moments technique. *Biomedical Engineering IEEE*, vol.50, 224–233.

Schwartz, M. Green, S. and Rutledge, W. 1960. *Vector Analysis with Applications to Geometry and Physics*. New York, USA: Harper and Brothers. 353p.

Thompson, L. and Pinsky, P. 1994. Complex wavenumber Fourier analysis of the p-version finite element method. *Computat. Mech.*, vol.13, 255–275.

Whitney, H. 1957. *Geometry Integration Theory*. New Jersey, USA: Princeton University Press. 341p.