I think I may try to formulate the problem using Lagrange interpolation to start off. I changed to a vector lagrange interpolation and MPC seemed to work. Question: I will need to invoke a Dirichlet condition on the perfect electric conducting surfaces. What is the best way to do \mathbf{n}\times\mathbf{E}=0 along a surface (where n is the normal unit vector along the Dirichlet boundary)? The boundaries are not generally aligned with the vector component directions.
UPDATE: I will try to use a weak form of the Dirichlet BCs, which I use in this version of the code.