HomerReid / scuff-em

A comprehensive and full-featured computational physics suite for boundary-element analysis of electromagnetic scattering, fluctuation-induced phenomena (Casimir forces and radiative heat transfer), nanophotonics, RF device engineering, electrostatics, and more. Includes a core library with C++ and python APIs as well as many command-line applications.
http://www.homerreid.com/scuff-em
GNU General Public License v2.0
126 stars 50 forks source link

Singularity Extraction Formulas #114

Closed urp1 closed 8 years ago

urp1 commented 8 years ago

Thanks of uploading this very useful code! I am having an issue using the standalone singularity program that I downloaded from (http://homerreid.dyndns.org/scuff-EM/SingularIntegrals/).

I am integrating 1/R over triangles T and T' using the code for T = T'. The results I obtain do not agree with the analytic solution - See the first entry in Table 1 in "On the Evaluation of the Double Surface Integrals Arising in the Application of the Boundary Integral Method to 3-D Problems, " Arcioni, Bressan, Perregrin, 1997 IEEE Transactions on MTT.

Scuff-EM works very well for all the cases I have tried so far, so I am not sure what is the reason for the mismatch. Test parameters I use with the code and the output are given below along with expected analytic solution.

I have set : int WhichCase = 3 double V1[3] = {0.0, 0.0, 0.0}; double V2[3] = {0.1, 0.0, 0.0}; double V3[3] = {0.05, 0.1, 0.0}; int NumPKs = 1; PIndex[1] = {TD_UNITY}; PIndex[1] = {TD_RP}; cdouble KParam[1] = {-1.0};

Result with code: 8.12533971e-1 Result with Analytic Formula: 1.0210603e-3

HomerReid commented 8 years ago

Thanks for your interest. The conversion factor here is 4*pi*(4*A*A') where A, A' are the triangle areas. In this case the areas are A=A'=0.005 so the conversion factor is 4*pi*4*0.005*0.005 ~= 0.00125, which explains your findings.

In other words, when you combine TD_UNITY with TD_RP, the actual integral you are getting is

\int dx \int dx' r^p / (4*pi*4*A*A')

I just checked and realized the online documentation is not quite right---the TD_RP kernel should read r^p/(4\pi) instead of r^p, while the TD_UNITY polynomial should read 1/(4AA^\prime) instead of 1. The former convention is designed in anticipation of using r^p integrals to desingularize the Helmholtz kernel, which has the 1/4\pi factor built in. The latter convention anticipates the use of RWG basis functions, for which the 1/(4AA^\prime) factor is always present by definition. I will update the documentation.

urp1 commented 8 years ago

Thanks for the quick response and fix!

urp1 commented 8 years ago

I noticed another issue with the singularity extraction routine when calculating the vector potential entries of the EFIE. Here are the parameters I have set:

int WhichCase = 3 double V1[3] = {0.0, 0.0, 0.0}; double V2[3] = {0.1, 0.0, 0.0}; double V3[3] = {0.05, 0.1, 0.0};

double Q[3] = {0.0, 0.0, 0.0}; double QP[3] = {0.01, 0.0, 0.0};

int NumPKs = 1; PIndex[1] = {TD_PMCHWG1}; PIndex[1] = {TD_RP}; cdouble KParam[1] = {-1.0};

"Scaled" result with the code: 2.4263e-7 Analytic result: -1.8635e-6 from formulas based on Arcioni, Bressan, and Perregrini paper and from Wilton, Rao, Glisson, Schaubert,Al-Bunduk, Butlers 1984 Transactions on AP paper.

When I use Q = QP = {0.0, 0.0, 0.0} the code works fine and I obtain the correct result of 4.07161e-6 with both "TaylorDuffy" routine and with the analytic formulas.

Is this a possible bug or I am making a mistake with the script?

Thanks!

HomerReid commented 8 years ago

The lines

PIndex[1] = {TD_PMCHWG1};
PIndex[1] = {TD_RP};

should read

PIndex[0] = {TD_PMCHWG1};
KIndex[0] = {TD_RP};

If you find a discrepancy with a result in the literature, please include a link to the reference so I can take a look.

Also, since the original issue you raised was solved, let's mark this issue as closed---if you find other problems or have other questions feel free to create new issues.

urp1 commented 8 years ago

I did use "KIndex". I have posted this as a different issue.

Thanks!!