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 on this Finite Difference Time Domain FDTD.  In this method, the space is discretized into a a grid of points and electromagnetic fields are evolved through time.  The next moment in time is calculated from the previous moment using only the information from the nearest points in space.  It is a direct application of 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.  In fact, this method will often not converge when non-physical items are introduced into the scene.

In contrast, the finite difference frequency domain (FD) techniques look at the scene from the point of view of a single frequency.  The most common variation on this is 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.  Using the frequency domain version of Maxwell's Equations, we can mathematically express how each point in space interacts with the other.  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 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.

(a) Permittivity as a function of space.  Black represents a permittivity of 4.  (b)  A plot of a snapshot of the phase of the electrical field.  (c)  The electric energy density.