Behaviour of DiffusionIntegrator when used as boundary integrator

[Rickbude]

Hi there, I have a question about the usage of a DiffusionIntegrator as a Boundary Integrator, and the behavior I should expect when doing so.

To understand why I even want to do this, let me describe the case first.

The case

I am trying to model "free-space"/"open" boundary conditions for the Poisson equation using Asymptotic Boundary conditions (ABC's), as explained in Gratkowski (2016) - "General Closed-Form Asymptotic Boundary Conditions for Finite Element Analysis of Exterior Electrical Field Problems" (doi: 10.15199/48.2016.05.04).

The weak form for the Poisson equation is given, after some manipulation, by:

$$\iint_{\Omega}{\epsilon \nabla T \cdot \nabla V dS} = \int_{\partial \Omega}{\epsilon T \nabla V \cdot \hat{n} dl} + \iint_{\Omega}{T\rho dS},$$

with $\epsilon$ the permittivity (let's take $\epsilon =\epsilon_0$ for simplicity), $V$ the potential, $T$ the test function. As an example, the second-order ABC boundary condition, expressed in 2D polar coordinates $(r,\phi)$, on the surface of a circle with $r=\mathrm{constant}=R$, is given by:

$$\nabla V \cdot \hat{n} = \frac{\partial V}{\partial r} = -\frac{1}{3R}\left(2V-\frac{\partial^2 V}{\partial \phi^2}\right)$$

This equation is obtained based on a multipole expansion of the potential. Substitute this into the weak form, and use partial integration to obtain for the boundary term (I think, I am not 100% sure about the partial integration):

$$\int_{\partial \Omega}{\epsilon T \nabla V \cdot \hat{n} dl} = - \frac{2\epsilon_0}{3R} \int_{\partial \Omega}{ TV dl} - \frac{\epsilon_0}{3R}\int_{\partial \Omega}{ \frac{\partial{T}}{\partial{\phi}}\frac{\partial{V}}{\partial{\phi}} dl}$$

Implementation options

Clearly the first term on the right-hand side can be implemented using a MassIntegrator (Or should I use a BoundaryMassIntegrator, is there a difference?). My confusion is mainly about the second term on the right-hand side. I can come up with two ways of implementing this:

Option 1

use $dl = rd\phi$ and the chain rule:

$$\frac{\partial{V}}{\partial\phi} = \frac{\partial{V}}{\partial{l}}\frac{\partial{l}}{\partial{\phi}}=r\frac{\partial{V}}{\partial{l}}$$

to write:

$$\int_{\partial \Omega}{\epsilon T \nabla V \cdot \hat{n} dl} = - \frac{2\epsilon_0}{3R} \int_{\partial \Omega}{ TV dl} - \frac{\epsilon_0R}{3}\int_{\partial \Omega}{ \frac{\partial{T}}{\partial{l}}\frac{\partial{V}}{\partial{l}} dl}$$

This seems to imply that I can use a DiffusionIntegrator with a scalar coefficient (or scaled identity matrix coefficient).

Option 2

Use the chain rule to write $\partial{V}/{\partial\phi}$ in terms of $\partial{V}/{\partial{x}}$ and $\partial{V}/{\partial{y}}$:

$$\frac{\partial{V}}{\partial\phi} = \frac{\partial{V}}{\partial{x}}\frac{\partial{x}}{\partial\phi} + \frac{\partial{V}}{\partial{y}}\frac{\partial{y}}{\partial\phi} = -y\frac{\partial{V}}{\partial{x}} +x\frac{\partial{V}}{\partial{y}}$$

and use that to write:

$$\int_{\partial \Omega}{\epsilon T \nabla V \cdot \hat{n} dl} = - \frac{2\epsilon_0}{3R} \int_{\partial \Omega}{ TV dl} - \frac{\epsilon_0R}{3}\int_{\partial \Omega}{ \left( \left[\begin{matrix} y^2 & -xy \\\ -xy & y^2 \end{matrix}\right] \nabla{V}\right) \cdot \nabla{T} dl }$$

Once again, we can use a DiffusionIntegrator, but now with a matrix coefficient. In contrast to the previous approach, this expands trivially to 3D simulations with 2D boundary segments.

Questions

While this type of boundary condition seems to work well or circular boundaries, the ellipsoidal and rectangular boundaries (for which a MixedScalarDerivativeIntegrator or MixedDirectionalDerivativeIntegrator is needed) have an error of 1-2 orders of magnitude higher on a comparable mesh size/density. I am trying to determine whether this is because I made an implementation or derivation error, or because I am maybe misunderstanding how bilinearform integrators involving the gradient work on the boundaries. This finally brings me to the questions:

• Does the DiffusionIntegrator, when used as a boundary integrator, operate in the same dimension number as the domain integrators? (2D in this case). Or on the dimension number of the boundary segment? (1D in this case). If the latter is the case, why is a 2x2 matrix coefficient accepted?
• Is option 1 or option 2 closer to how the code works? Which option is "better"? Strangely enough, the results are exactly the same in this example case (See picture below).

My apologies for this long post and the many questions. Thanks in advance!