Finite Element Method for Time-Harmonic Maxwell’s Equations#

This section presents two finite-element formulations of the time-harmonic Maxwell equations. We first introduce the conventional weak form (the primal formulation), then the ultraweak formulation used by Lethe’s DPG solver.

Primal Formulation#

Starting from the electric-field wave equation shown in Time Harmonic Maxwell’s Equations:

\[\nabla \times \left( \mu_{\mathrm{em}}^{-1} \nabla \times \mathbf{E} \right) -\omega^2 \varepsilon_{\mathrm{em,eff}} \mathbf{E} = i \omega \mathbf{J}_{\mathrm{ext}},\]

we consider a domain \(\Omega\) with boundary \(\Gamma\).

Note

In general, it is not necessary to solve both the electric- and magnetic-field wave equations, since one can be derived from the other. The field to solve for can be chosen based on the problem at hand.

Boundary conditions depend on the problem (for example Sommerfeld, absorbing/radiation, or impedance conditions), but they can be grouped into three main categories:

First-type (Dirichlet) boundary conditions: \(\mathbf{n} \times \mathbf{E} = \mathbf{E}_{\mathrm{D}}\);
Second-type (Neumann) boundary conditions: \(\mu_{\mathrm{em}}^{-1} \mathbf{n} \times (\nabla \times \mathbf{E}) = i \omega \mathbf{J}_{\mathrm{ext,N}}\);
Robin boundary conditions: \(\mathbf{n} \times (\mu_{\mathrm{em}}^{-1} \nabla \times \mathbf{E}) + Y_\mathrm{s} \mathbf{n} \times ( \mathbf{E} \times \mathbf{n}) = \mathbf{E}_\mathrm{R}\),

In the expressions above, \(\mathbf{n}\) is the outward unit normal and \(Y_\mathrm{s}\) is the boundary surface admittance. The Robin condition generalizes Dirichlet and Neumann conditions, and can represent impedance, absorbing, and related boundary models.

For simplicity, we consider a perfect electric conductor (PEC), so \(\mathbf{E}_{\mathrm{D}} = 0\). Multiplying the strong form by a complex test function \(\mathbf{v}\) satisfying \(\mathbf{v} \times \mathbf{n}=0\), and integrating over \(\Omega\), gives:

\[\begin{split}\begin{align*} \int_{\Omega} i \omega \mathbf{J}_{\mathrm{ext}} \cdot \mathbf{v^*} \mathrm{d}\Omega = & \int_{\Omega} \mu_{\mathrm{em}}^{-1} (\nabla \times \mathbf{E}) \cdot (\nabla \times \mathbf{v^*}) \mathrm{d}\Omega - \int_{\Omega} \omega^2 \varepsilon_{\mathrm{em},\mathrm{eff}} \mathbf{E} \cdot \mathbf{v^*} \mathrm{d}\Omega \\ & + \int_{\Gamma} \mu_{\mathrm{em}}^{-1} (\nabla \times \mathbf{E}) \cdot (\mathbf{v^*} \times \mathbf{n}) \mathrm{d}\Gamma . \end{align*}\end{split}\]

Because the fields are complex-valued, the usual \(L^2\) inner product is replaced by \(\int_{\Omega} \mathbf{u}\mathbf{v^*} \, \mathrm{d}\Omega\), where \(\mathbf{v^*}\) denotes the complex conjugate of \(\mathbf{v}\). For these integrals to be well-defined, \(\mathbf{E}\) and \(\mathbf{v}\) must belong to the Sobolev space \(H(\mathrm{curl}, \Omega)\):

\[H(\mathrm{curl}, \Omega) = \{ \mathbf{v} \in [L^2(\Omega)]^3 : \nabla \times \mathbf{v} \in [L^2(\Omega)]^3 \},\]

and the Dirichlet boundary condition imposes the homogeneous curl space:

\[H_0(\mathrm{curl}, \Omega) = \{ \mathbf{v} \in H(\mathrm{curl}, \Omega) : \mathbf{n} \times \mathbf{v}|_{\Gamma} = 0 \}.\]

Thus, the boundary term vanishes (\(\int_{\Gamma} \mu_{\mathrm{em}}^{-1} (\nabla \times \mathbf{E}) \cdot (\mathbf{v^*} \times \mathbf{n}) \mathrm{d}\Gamma = 0\)), and the weak form of the time-harmonic Maxwell equations becomes:

\[\begin{split}B(\mathbf{E}, \mathbf{v}) = &\int_{\Omega} \mu_{\mathrm{em}}^{-1} (\nabla \times \mathbf{E}) \cdot (\nabla \times \mathbf{v^*}) \mathrm{d}\Omega - \int_{\Omega} \omega^2 \varepsilon_{\mathrm{em},\mathrm{eff}} \mathbf{E} \cdot \mathbf{v^*} \mathrm{d}\Omega \\ L(\mathbf{v}) = &\int_{\Omega} i \omega \mathbf{J}_{\mathrm{ext}} \cdot \mathbf{v^*} \mathrm{d}\Omega .\end{split}\]

Formally, \(\mathbf{E}\) should also satisfy Gauss’s law (\(\nabla \cdot \mathbf{D} = \rho_f\)). In this setting, it is implicitly enforced by the electromagnetic wave equation and is therefore contained in the weak form above.

Ultraweak Formulation#

The formulation above is the conventional weak form, used here as a reference for comparison with other implementations of the time-harmonic Maxwell equations. In Lethe’s DPG setting, the ultraweak formulation is used instead. The main practical difference is that the primal formulation above eliminates the magnetic field to obtain a single second-order equation for the electric field, whereas the ultraweak formulation solves for the electric and magnetic fields simultaneously by weakening both equations without eliminating either one.

Retaining both fields as independent unknowns is especially useful when the electromagnetic power dissipated in the medium must be evaluated, since some materials exhibit non-negligible magnetic losses in addition to dielectric losses; recovering \(\mathbf{H}\) directly, rather than reconstructing it from \(\mathbf{E}\) after the fact, avoids an additional loss of accuracy in that calculation. Additionally, studies have shown that the ultraweak formulation is more robust than the primal formulation in the context of the DPG method.

The starting point is the time-harmonic Maxwell system written in the form

\[\begin{split}\begin{align} \nabla \times \mathbf{E} - i\omega \mu_{r} \mathbf{H} &= 0, \\ \nabla \times \mathbf{H} + i\omega \varepsilon_{r,\mathrm{eff}} \mathbf{E} &= \mathbf{J}, \end{align}\end{split}\]

where \(\varepsilon_{r,\mathrm{eff}} = \varepsilon_r + i \frac{\sigma_e}{\omega \varepsilon_0}\). The system is supplemented with boundary conditions of three kinds:

\[\begin{split}\begin{align} \mathbf{n} \times \mathbf{E} &= \mathbf{n} \times \mathbf{E}_D; \\ \mathbf{n} \times \mathbf{H} &= \mathbf{n} \times \mathbf{H}_N - \mathbf{J}_{\mathrm{s}}; \\ \mathbf{n} \times \mathbf{H} + Z_{\mathrm{s}}^{-1}\mathbf{n} \times (\mathbf{E}\times \mathbf{n}) &= \mathbf{E}_\mathrm{R}. \end{align}\end{split}\]

To obtain the ultraweak formulation, Faraday’s law and Ampère-Maxwell equations are multiplied by independent test functions \(\mathbf{I}\) and \(\mathbf{F}\), respectively, and integrated by parts over the mesh \(\Omega_h\):

\[\begin{split}\begin{align} (\nabla \times\mathbf{I},\mathbf{E})_{\Omega_h} - ( \mathbf{I}, i\omega \mu_r\mathbf{H})_{\Omega_h} + \langle \mathbf{I}, \mathbf{n}\times \hat{\mathbf{E}} \rangle_{\partial \Omega_h} &= 0 , \\ (\nabla \times \mathbf{F},\mathbf{H})_{\Omega_h} + (\mathbf{F}, i \omega \varepsilon_{r,\text{eff}} \mathbf{E})_{\Omega_h} + \langle \mathbf{F} , \mathbf{n}\times \hat{\mathbf{H}}\rangle_{\partial \Omega_h \backslash \Gamma_R } - \langle \mathbf{F} ,Y_\mathrm{s} \mathbf{n} \times ( \mathbf{E} \times \mathbf{n})\rangle_{ \Gamma_R } &=(\mathbf{F},\mathbf{J})_{\Omega_h} - \langle \mathbf{F} , \mathbf{E_\mathrm{R}} \rangle_{ \Gamma_R }, \end{align}\end{split}\]

where the Robin boundary condition is applied to the electric field without loss of generality, and \(\hat{\mathbf{E}}\), \(\hat{\mathbf{H}}\) are the trace unknowns on the mesh skeleton. Each product above is understood to be complex-conjugate in its first argument:

\[\begin{align} (\mathbf{v},\mathbf{u})_{\Omega_h} &:= \sum_{K\in\Omega_h} \int_{K} \mathbf{v}^* \cdot \mathbf{u}\,\mathrm{d}x, & \langle \mathbf{v},\mathbf{u} \rangle_{\partial\Omega_h} &:= \sum_{K\in\Omega_h} \int_{\partial K} \mathbf{v}^* \cdot \mathbf{u}\, \mathrm{d}s, \end{align}\]

with \(\mathbf{v}\) and \(\mathbf{u}\) being arbitrary complex-valued vector fields. From the above ultraweak system, the trial and test spaces are defined as :

\[\begin{split}\begin{align} \mathbf{E},\mathbf{H} &\in \mathcal{X}_{\mathrm{em}}, & \mathcal{X}_{\mathrm{em}} &:= (L^2(\Omega))^3, \\ \hat{\mathbf{E}},\hat{\mathbf{H}} &\in \hat{\mathcal{X}}_{\mathrm{em}}, & \hat{\mathcal{X}}_{\mathrm{em}} &:= \bigl\{ \hat{\mathbf{x}} \in H^{-1/2}(\mathrm{curl},\partial\Omega_h) : \mathbf{n}\times \hat{\mathbf{x}} = \mathbf{n}\times \mathbf{x}_\mathrm{D} \text{ on } \Gamma_D \bigr\}, \\ \mathbf{F},\mathbf{I} &\in \mathcal{Y}_{\mathrm{em}}, & \mathcal{Y}_{\mathrm{em}} &:= \bigl\{ \mathbf{y} \in H(\mathrm{curl},\Omega_h) \bigr\}. \end{align}\end{split}\]

For an ultraweak formulation, the Neumann boundary condition is, in fact, a Dirichlet condition on the flux variable (here, the magnetic field \(\mathbf{H}\)). This is why only a Dirichlet-type condition needs to be imposed on the trace space \(\hat{\mathcal{X}}_{\mathrm{em}}\), which covers both the electric and magnetic trace unknowns.

To complete the discrete DPG formulation, the test space is equipped with the adjoint graph norm induced by the operator associated with the bilinear form of the system:

\[\begin{split}\begin{align} \|(\mathbf{F}, \mathbf{I})\|^2_V &= \| \nabla \times \mathbf{F} - i \omega \mu_r \mathbf{I} \|^2_{\Omega_h} + \| \nabla \times \mathbf{I} + i \omega \varepsilon_{r,\mathrm{eff}} \mathbf{F} \|^2_{\Omega_h} \\ &\quad + \alpha\left(\|\mathbf{F}\|^2_{\Omega_h} + \|\mathbf{I}\|^2_{\Omega_h}\right) + \| \mathbf{n} \times \mathbf{I} + Y_\mathrm{s} \mathbf{n} \times (\mathbf{F} \times \mathbf{n}) \|^2_{\Gamma_R}, \end{align}\end{split}\]

where \(\alpha \sim 1\) is a strictly positive constant (\(\alpha=1\) in practice) introduced to ensure that the test norm remains localizable once the test space is broken. The boundary term on \(\Gamma_R\) enforces the Robin boundary condition and is required to obtain convergence of the method when this type of boundary condition is present.

Tip

The general Dirichlet, Neumann, and Robin conditions introduced above cover a range of standard boundary condition types encountered in electromagnetic simulations, each corresponding to a particular choice of \(\mathbf{E}_\mathrm{D}\), \(\mathbf{H}_\mathrm{N}\), \(Y_\mathrm{s}\), and \(\mathbf{E}_\mathrm{R}\). Writing \(\mathbf{E}_\parallel = \mathbf{n} \times (\mathbf{E} \times \mathbf{n})\) for the tangential component of a field (and similarly \(\mathbf{H}_\parallel\)), the list below is not exhaustive but covers the conditions commonly used in practice.

Perfect electric conductor (PEC): \(\mathbf{n} \times \mathbf{E} = 0\).
Perfect magnetic conductor (PMC): \(\mathbf{n} \times \mathbf{H} = 0\).
Prescribed electric field (Dirichlet): \(\mathbf{n} \times \mathbf{E} = \mathbf{n} \times \mathbf{E}_D\).
Prescribed magnetic field (Neumann): \(\mathbf{n} \times \mathbf{H} = \mathbf{n} \times \mathbf{H}_N - \mathbf{J}_{\mathrm{s}}\), with \(\mathbf{J}_{\mathrm{s}}\) a prescribed surface current density.
Impedance condition (Robin): \(\mathbf{n} \times \mathbf{H} + Y_\mathrm{s} \mathbf{E}_\parallel = \mathbf{E_\mathrm{R}}\), with \(Y_\mathrm{s}\) the surface admittance of the boundary.
Absorbing / radiation condition (Silver–Müller): \(\mathbf{n} \times \mathbf{H} + \sqrt{\varepsilon_{r,\mathrm{eff}}/\mu_r}\, \mathbf{E}_\parallel = 0\). This is the impedance condition above with \(Y_\mathrm{s} = \sqrt{\varepsilon_{r,\mathrm{eff}}/\mu_r}\) (the admittance of the exterior medium) and no source term.
Lossy / imperfect conductor: \(\mathbf{n} \times \mathbf{H} + \sqrt{\varepsilon_{r,\mathrm{eff},2}/\mu_{r,2}}\, \mathbf{E}_\parallel = 0\), the same form as the radiation condition above, using the (possibly complex) effective properties \(\varepsilon_{r,\mathrm{eff},2}\), \(\mu_{r,2}\) of the adjacent conducting medium.
Waveguide inlet port, TE\(_{mn}\) mode: \(\mathbf{n} \times \mathbf{H} + \dfrac{\mathbf{k} \cdot \mathbf{n}}{\omega \mu_r} \mathbf{E}_\parallel = \mathbf{n} \times \mathbf{H}_{inc} + \dfrac{\mathbf{k} \cdot \mathbf{n}}{\omega \mu_r} \mathbf{E}_{inc,\parallel}\).
Waveguide inlet port, TM\(_{mn}\) mode: \(\mathbf{n} \times \mathbf{H} + \dfrac{\omega \varepsilon_{r,\mathrm{eff}}}{\mathbf{k} \cdot \mathbf{n}} \mathbf{E}_\parallel = \mathbf{n} \times \mathbf{H}_{inc} + \dfrac{\omega \varepsilon_{r,\mathrm{eff}}}{\mathbf{k} \cdot \mathbf{n}} \mathbf{E}_{inc,\parallel}\).

In each waveguide, \(\mathbf{H}_{inc}\) and \(\mathbf{E}_{inc,\parallel}\) denote the incident field of the excited mode, and \(k = \omega \sqrt{\varepsilon_{r,\mathrm{eff}} \mu_r}\) its wavevector amplitude.