Finite Difference Frequency Domain in Mathematica

Electromagnetic problems are typically solved using one of two classes of numerical solvers.  The first is time-domain (TD).  The most common variation of this Finite Difference Time Domain (FDTD).  In this method, space is discretized into a grid of points, and electromagnetic fields evolve over time.  The next moment in time is calculated from the previous moment using only the information from the nearest points in space.  It directly applies the time-domain of Maxwell's equations and reflects principles of locality in both space and causality in time.  It is very similar to how nature "calculates" electromagnetic fields.  This method will often not converge when non-physical items are introduced into the scene.

In contrast, the frequency domain (FD) techniques look at the scene from the point of view of a single frequency.  The most common variation on this is the Finite Element Method (FEM).  FEM discretizes the scene into volumetric elements using polynomial functions and then seeks to minimize the unresolved energy.  The Finite Difference Frequency Domain (FDFD) technique uses a similar hexahedral mesh as the FDTD technique, but also the frequency domain version of Maxwell's Equations. This formulation expresses how each point in space interacts with its neiboring points.  This creates a system of constraints which, if formulated correctly, will be equal in size to the number of unknowns in the system.  In other words, the linear system is determined and square.  Solving this problem reduces to inverting this system.

I implemented the FDFD technique in Mathematica.  The permittivity can be specified as a function of geometry.  Additionally, it supports PEC, PMC and PML.  An example output is below.

I have also implemented the FDFD technique in Python.