Singularity Routine

[urp1]

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 main(int argc, char *argv[])
{
  /* panel vertices */
  int WhichCase = 3; // 2 common vertices
  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.1, 0.0, 0.0};

  /* specification of which integrals we want */
  int NumPKs        = 1;
  int PIndex[1]     = {TD_PMCHWG1};
  int KIndex[1]     = {TD_RP};
  cdouble KParam[1] = {-1.0};

  /* output buffers */
  cdouble Result[1], Error[1];

  /* fill in argument structure */
  TaylorDuffyArgStruct MyTDArgs, *TDArgs=&MyTDArgs;
  InitTaylorDuffyArgs(TDArgs);
  TDArgs->WhichCase = WhichCase;
  TDArgs->V1        = V1;
  TDArgs->V2        = V2;
  TDArgs->V3        = V3;
  TDArgs->Q         = Q;
  TDArgs->QP        = QP;
  TDArgs->NumPKs    = NumPKs;
  TDArgs->PIndex    = PIndex;
  TDArgs->KIndex    = KIndex;
  TDArgs->KParam    = KParam;
  TDArgs->Result    = Result;
  TDArgs->Error     = Error;  

  /* evaluate the integrals */
  TaylorDuffy( TDArgs );

  double A = 0.005; //Area of unprimed triangle
  double Ap = 0.005;    //Area of primed triangle
  double PI = 3.14159265359;

  cdouble Result_1 = Result[0]*(4.0*A*Ap)*4.0*PI;
  std::cout << Result_1 << std::endl;

}

Executing the code: 2.4263e-7

Analytic result: -1.8635e-6 obtained with the formulas in the following paper: http://ieeexplore.ieee.org/document/563344/?arnumber=563344&tag=1

Here is the c++ implementation of the code for the formulas in this paper

double SingularityAnalytic2(double a, double b, double c, double p, double A, bool flag)
{
    // a: length 1
    // b: length 2
    // c: length 3
    // p: half-perimeter
    // A: Area
    // flag: 1 if source and test basis are the same.
    // flag: 0 if source and test basis are different.

    double Result;
    double a2 = pow(a,2);
    double b2 = pow(b,2);
    double c2 = pow(c,2);
    double A2 = pow(A,2);
    if(flag)
    {
        Result = A2/(30.0)*((10.0 + 3.0*(c2-a2)/b2 - 3.0*(a2-b2)/c2)*a - (5.0-3.0*(a2-b2)/c2 - 2.0*(b2-c2)/a2)*b - (5.0+3.0*(c2-a2)/b2 + 2.0*(b2-c2)/a2)*c
            + (a2 - 3.0*b2 - 3.0*c2 - 8.0*A2/a2)*2.0/a*log(1-a/p) + (a2-2.0*b2-4.0*c2+6.0*A2/b2)*4.0/b*log(1.0-b/p) + (a2-4.0*b2 - 2.0*c2 + 6.0*A2/c2)*4.0/c*log(1-c/p));

    }
    else{
        //Case alpha not equal to beta
        Result = A2/60.0*((-10.0 + (c2-a2)/b2 - (a2-b2)/c2)*a + (5.0 + (a2-b2)/c2 - 6.0*(b2-c2)/a2)*b + (5.0 - (c2-a2)/b2 + 6.0*(b2-c2)/a2)*c
                + (2.0*a2 - b2 - c2 + 4.0*A2/a2)*12.0/a*log(1.0-a/p) + (9.0*a2 - 3.0*b2 - c2 + 4.0*A2/b2)*2.0/b*log(1.0-b/p) + (9.0*a2 - b2 - 3.0*c2 + 4.0*A2/c2)*2.0/c*log(1.0-c/p));
    }

    return Result;
}

Which is executed by calling the following:

    double a = 0.1;
    double b = sqrt(0.1*0.1 + 0.05*0.05);
    double c = sqrt(0.1*0.1 + 0.05*0.05);
    double p = (a+b+c)/2;
    double A = 0.005;
    double Result1 = SingularityAnalytic2(a, b, c, p, A, true);         //Basis FV is same as testing FV.

    double Result2 = SingularityAnalytic2(a,b,c, p, A, false);

    cout << "Result 1 " << Result1 << endl;
    cout << "Result 2 " << Result2 << endl;

The singularity routine with scuff-em does provide correct result compared to analytic formulas when alpha = beta, i.e. when the Q = QP.

[HomerReid]

You've incorrectly identified QP as coinciding with vertex V2 in your call to my code, whereas in fact it seems QP coincides with V3 in the example against which you are comparing.

If I modify your call to my code to read

double QP[3] = {0.05, 0.1, 0.0};

then I get -1.86352e-6 which is the answer you want.

[urp1]

Yes you are right. The result with Analytic formula for that case was -1.03e-6, I jumbled the vertex.

Nevertheless, if set QP to the value you indicated then I get: (1.96546e-06,0).

Is it possible that the code on the website (http://homerreid.dyndns.org/scuff-EM/SingularIntegrals/) is different then the code that you have? Perhaps some issues were resolved afterwards?

I appreciate your help.

[HomerReid]

I don't think so, but you can easily test this by downloading and installing SCUFF-EM and then compiling the test code against libscuff. To do this, add the following two lines to the top of the test code:

#include <TaylorDuffy.h>
using namespace scuff;

Then compile and link with flags such as

CPPFLAGS=-I/path/to/scuff-em-installation/include/scuff-em
LDFLAGS=-L/path/to/scuff-em-installation/lib

with -lscuff added to the linker command line.

[urp1]

The version in scuff-em package gives correct results! So there must be some bug in the Singularity extraction package on the website.

Thanks for your help!

[HomerReid]

Thanks for identifying this issue.

I updated the website documentation to eliminate the old buggy code and to clean up the notational confusion that was the subject of your original issue.

Let me know if you have any further questions.

Also, if you are interested, I have an upcoming publication presenting a Taylor-Duffy method for tetrahedron-product integrals (for volume-integral-equation solvers).

[urp1]

Thanks!

Yes, I am interested in reading your upcoming publication, where can I download it? Singularity handling is always tricky, but your package is very useful.

[HomerReid]

The publication will be up on ArXiV next week. The code is already available for download in my BUFF-EM package but is not as thoroughly documented as in SCUFF-EM.