Working on the electromagnetic version of the time-harminic solver. I am using the first order Maxwell curl equations in an impedance-normalized form (i.e. \eta_0 = \sqrt{\mu_0 / \epsilon_0} = 1). I consider the case of inhomogeneous dielectrics only (hence \epsilon_r in the following equations). The wavenumber (spatial frequency) is k_0.
\nabla\times\mathbf{E} = -jk_0 \mathbf{H}
\nabla\times\mathbf{H} = jk_0\epsilon_r\mathbf{E}
Like in the acoustic case, I form the least-squares residual using the canonical first-order system and the relevant boundary conditions:
R = \vert\!\vert\nabla\times\mathbf{E}+jk_0\mathbf{H}\vert\!\vert^2_\Omega + \vert\!\vert\nabla\times\mathbf{H}-jk_0\epsilon_r\mathbf{E}\vert\!\vert^2_\Omega + \vert\!\vert\mathbf{n\times E}\vert\!\vert^2_{\partial\Omega_{PEC}}+\vert\!\vert\mathbf{n\cdot H}\vert\!\vert^2_{\partial\Omega_{PEC}} +
\hskip{1in}\vert\!\vert\mathbf{n\times ((E - E^{inc})\times n)} - Z_b \mathbf{ n\times H)}\vert\!\vert^2_{\partial\Omega_{InOut}}
Z_b is a problem-dependent normalized (to \eta_0) boundary impedance used to absorb outgoing waves. I am using curl-conforming Nédélec functions as the finite element interpolation space. Preliminary results look good, but I need to add more memory to my laptop to verify convergence! 
The boundary conditions are applied quite differently than in the usual Galerkin method. Moreover, low order interpolation tends to exhibit a non-physical loss. At least order 2 interpolation is required for decent performance. Order 3 is best, but generates huge matrix orders (with higher density as well).
Until next time!