add model HelmholtzHABCwithCorners

HelmholtzHABCwithCorners/README.txt

+High-order absorbing boundary conditions for 2D and 3D Helmholtz problems.
+
+A. Modave, X. Geuzaine, X. Antoine (2020). Corner treatments for high-order absorbing boundary conditions in high-frequency acoustic scattering problems. J. Comput. Phys., 401, 109029 (preprint: https://hal.archives-ouvertes.fr/hal-01925160)
+
+Models developed by Axel Modave.
+
+
+Quick start
+-----------
+
+Open `main.pro' with Gmsh.
\ No newline at end of file

HelmholtzHABCwithCorners/main.dat

+DIR = "out/";
+
+GEO_SQUARE     = 1;
+GEO_DISK       = 2;
+GEO_POLYGON    = 3;
+GEO_SLICE      = 4;
+GEO_CUBE       = 5;
+GEO_POLYHEDRON = 6;
+SIGNAL_Harmonic  = 1;
+SIGNAL_Scatt     = 2;
+SIGNAL_Dirichlet = 1;
+SIGNAL_Neumann   = 2;
+
+DefineConstant[
+  FLAG_GEO = {GEO_SQUARE,
+    Name "Input/1Geometry/0Type of domain", Highlight "Blue",
+    GmshOption "Reset", Autocheck 0,
+    Choices {GEO_SQUARE     = "Squared domain (2D)",
+             GEO_DISK       = "Circular domain (2D)",
+             GEO_POLYGON    = "Polygonal domain (2D)",
+             GEO_SLICE      = "Slice domain (2D)",
+             GEO_CUBE       = "Cuboidal domain (3D)",
+             GEO_POLYHEDRON = "Polyhedral domain (3D)"}}
+];
+
+If ((FLAG_GEO==GEO_SQUARE) || (FLAG_GEO==GEO_DISK) || (FLAG_GEO==GEO_POLYGON) || (FLAG_GEO==GEO_SLICE))
+  FLAG_DIM = 2;
+  defWaveNum = 25;
+  defNLambda = 10;
+ElseIf ((FLAG_GEO==GEO_CUBE) || (FLAG_GEO==GEO_POLYHEDRON))
+  FLAG_DIM = 3;
+  defWaveNum = 10;
+  defNLambda = 10;
+EndIf
+
+DefineConstant[
+  FLAG_SIGNAL = {SIGNAL_Scatt,
+    Name "Input/2Signal/1Type of solution",
+    Choices {SIGNAL_Harmonic = "Harmonic (cylindrical or spherical)",
+             SIGNAL_Scatt    = "Scattering of plane wave"}},
+  FLAG_SIGNAL_BC = {SIGNAL_Neumann,
+    Name "Input/2Signal/2Type of condition",
+    Choices {SIGNAL_Dirichlet = "Sound-soft (Dirichlet BC)",
+             SIGNAL_Neumann   = "Sound-hard (Neumann BC)"}},
+  SIGNAL_MODE = {0,          Min 0,   Step 1,   Max 50,              Name "Input/2Signal/3Mode number", Visible (FLAG_SIGNAL == SIGNAL_Harmonic)},
+  WAVENUMBER  = {defWaveNum, Min 0.1, Step 0.1, Max 300,             Name "Input/2Signal/4Wavenumber"},
+  LAMBDA      = {2*Pi/WAVENUMBER,                                    Name "Input/2Signal/5Wavelength", ReadOnly 1},
+  R_SCA       = {1,          Min 0.1, Step 0.1, Max  10,             Name "Input/1Geometry/8Scatterer radius"},
+  N_LAMBDA    = {defNLambda, Choices {5, 7.5, 10, 12.5, 15},         Name "Input/3Discretization/1Points per wavelength"},
+  ORDER       = {2, Choices {1 = "First-order", 2 = "Second-order"}, Name "Input/3Discretization/2Polynomial order"}
+];
+
+// LAMBDA = 2.*Pi/WAVENUMBER;
+// WAVENUMBER = 2.*Pi/LAMBDA;
+LC = LAMBDA/N_LAMBDA;
+
+LinkGeo = 0;
+LinkPro = 0;
+If (FLAG_GEO==GEO_SQUARE)
+  LinkGeo = "padeSquare.geo";
+  LinkPro = "padeSquare.pro";
+ElseIf (FLAG_GEO==GEO_DISK)
+  LinkGeo = "padeDisk.geo";
+  LinkPro = "padeDisk.pro";
+ElseIf (FLAG_GEO==GEO_POLYGON)
+  LinkGeo = "padePolygon.geo";
+  LinkPro = "padePolygon.pro";
+ElseIf (FLAG_GEO==GEO_SLICE)
+  LinkGeo = "padePieWedge.geo";
+  LinkPro = "padePieWedge.pro";
+ElseIf (FLAG_GEO==GEO_CUBE)
+  LinkGeo = "padeCube.geo";
+  LinkPro = "padeCube.pro";
+ElseIf (FLAG_GEO==GEO_POLYHEDRON)
+  LinkGeo = "padePolyhedron.geo";
+  LinkPro = "padePolyhedron.pro";
+EndIf

HelmholtzHABCwithCorners/main.geo

+Include "main.dat" ;
+
+CreateDir Str(DIR);
+
+SetOrder ORDER;
+Mesh.ElementOrder = ORDER;
+Mesh.SecondOrderLinear = 0;
+
+Mesh.CharacteristicLengthMax = LC;
+Mesh.CharacteristicLengthFactor = 1;
+Mesh.Optimize = 1;
+
+// Solver.AutoMesh = 1;
+
+Include Str[LinkGeo];

HelmholtzHABCwithCorners/main.pro

+Include "main.dat" ;
+
+//==================================================================================================
+// FUNCTIONS FOR SIGNAL
+//==================================================================================================
+
+Function {
+  I[] = Complex[0,1];
+  k[] = WAVENUMBER;
+  R[] = Sqrt[X[]^2+Y[]^2+Z[]^2];
+  THETA[] = Atan2[Y[],X[]];
+
+If((FLAG_DIM == 2) && (FLAG_SIGNAL == SIGNAL_Harmonic))
+  ExpMode[] = Complex[Cos[SIGNAL_MODE*THETA[]], Sin[SIGNAL_MODE*THETA[]]];
+  JnR[] = Jn[SIGNAL_MODE,k[]*R_SCA];
+  Jnr[] = Jn[SIGNAL_MODE,k[]*R[]  ];
+  YnR[] = Yn[SIGNAL_MODE,k[]*R_SCA];
+  Ynr[] = Yn[SIGNAL_MODE,k[]*R[]  ];
+  dJnR[] = k[]*dJn[SIGNAL_MODE,k[]*R_SCA];
+  dJnr[] = k[]*dJn[SIGNAL_MODE,k[]*R[]  ];
+  dYnR[] = k[]*dYn[SIGNAL_MODE,k[]*R_SCA];
+  dYnr[] = k[]*dYn[SIGNAL_MODE,k[]*R[]  ];
+  HnkR[] = Complex[JnR[], YnR[]];
+  Hnkr[] = Complex[Jnr[], Ynr[]];
+  dHnkR[] = Complex[dJnR[], dYnR[]];
+  dHnkr[] = Complex[dJnr[], dYnr[]];
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+  f_inc[] = Jnr[] * ExpMode[];
+  f_ref[] = -(JnR[]/HnkR[]) * Hnkr[] * ExpMode[];
+EndIf
+If(FLAG_SIGNAL_BC == SIGNAL_Neumann)
+  df_inc[] = dJnr[] * ExpMode[];
+  f_ref[] = -(dJnR[]/dHnkR[]) * Hnkr[] * ExpMode[];
+EndIf
+EndIf
+
+If((FLAG_DIM == 2) && (FLAG_SIGNAL == SIGNAL_Scatt))
+  f_inc[] = Complex[Cos[k[]*X[]], Sin[k[]*X[]]];
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+  f_ref[] = AcousticFieldSoftCylinder[XYZ[]]{WAVENUMBER, R_SCA};
+EndIf
+If(FLAG_SIGNAL_BC == SIGNAL_Neumann)
+  df_inc[] = k[] * Complex[-Sin[k[]*X[]], Cos[k[]*X[]]] * (X[]/R[]);  // Grad * \hat{r}
+  // df2_inc[] = k[] * Complex[-Sin[k[]*X[]], Cos[k[]*X[]]];             //  CompX[Grad]
+  f_ref[] = AcousticFieldHardCylinder[XYZ[]]{WAVENUMBER, R_SCA, 0, 0, 50};
+EndIf
+EndIf
+
+If((FLAG_DIM == 3) && (FLAG_SIGNAL == SIGNAL_Harmonic))
+  JnSphR[] = JnSph[SIGNAL_MODE,k[]*R_SCA];
+  JnSphr[] = JnSph[SIGNAL_MODE,k[]*R[]  ];
+  YnSphR[] = YnSph[SIGNAL_MODE,k[]*R_SCA];
+  YnSphr[] = YnSph[SIGNAL_MODE,k[]*R[]  ];
+  dJnSphR[] = k[]*dJnSph[SIGNAL_MODE,k[]*R_SCA];
+  dJnSphr[] = k[]*dJnSph[SIGNAL_MODE,k[]*R[]  ];
+  dYnSphR[] = k[]*dYnSph[SIGNAL_MODE,k[]*R_SCA];
+  dYnSphr[] = k[]*dYnSph[SIGNAL_MODE,k[]*R[]  ];
+  HnSphkR[] = Complex[JnSphR[], YnSphR[]];
+  HnSphkr[] = Complex[JnSphr[], YnSphr[]];
+  dHnSphkR[] = Complex[dJnSphR[], dYnSphR[]];
+  dHnSphkr[] = Complex[dJnSphr[], dYnSphr[]];
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+  f_inc[] = JnSphr[];
+  f_ref[] = -(JnSphR[]/HnSphkR[]) * HnSphkr[];
+EndIf
+If(FLAG_SIGNAL_BC == SIGNAL_Neumann)
+  df_inc[] = dJnSphr[];
+  f_ref[] = -(dJnSphR[]/dHnSphkR[]) * HnSphkr[];
+EndIf
+EndIf
+
+If((FLAG_DIM == 3) && (FLAG_SIGNAL == SIGNAL_Scatt))
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+  f_inc[] = Complex[Cos[k[]*X[]], Sin[k[]*X[]]];
+  f_ref[] = AcousticFieldSoftSphere[XYZ[]]{WAVENUMBER, R_SCA, 1, 0, 0};
+EndIf
+If(FLAG_SIGNAL_BC == SIGNAL_Neumann)
+  df_inc[] = Complex[0,1] * k[]*X[]/R[] * Complex[Cos[k[]*X[]], Sin[k[]*X[]]];
+  f_ref[] = AcousticFieldHardSphere[XYZ[]]{WAVENUMBER, R_SCA, 1, 0, 0};
+EndIf
+EndIf
+}
+
+//==================================================================================================
+// JACOBIAN & INTEGRATION
+//==================================================================================================
+
+Jacobian {
+  { Name JVol; Case {{ Region All; Jacobian Vol; }}}
+  { Name JSur; Case {{ Region All; Jacobian Sur; }}}
+  { Name JLin; Case {{ Region All; Jacobian Lin; }}}
+}
+
+Integration {
+  { Name I1;
+    Case {
+      { Type Gauss;
+        Case {
+          { GeoElement Point;        NumberOfPoints  1; }
+          { GeoElement Line;         NumberOfPoints  6; } // 4
+          { GeoElement Line2;        NumberOfPoints  6; } // 4
+          { GeoElement Quadrangle;   NumberOfPoints  4; } // 36
+          { GeoElement Quadrangle2;  NumberOfPoints  7; } // 36
+          { GeoElement Triangle;     NumberOfPoints  4; } // 6 // 1 3 4 6 7 12 13 16
+          { GeoElement Triangle2;    NumberOfPoints  6; } // 6 // 1 3 4 6 7 12 13 16 // 12
+          { GeoElement Tetrahedron;  NumberOfPoints  5; } // 15 // 1 4 5 15 16 17 29
+          { GeoElement Tetrahedron2; NumberOfPoints 15; } // 15 // 1 4 5 15 16 17 29
+        }
+      }
+    }
+  }
+}
+
+//==================================================================================================
+// LOAD SPECIFIC .PRO
+//==================================================================================================
+
+Include Str[LinkPro];
+
+DefineConstant[
+  C_ = {"-solve -pos -bin -v2", Name "GetDP/9ComputeCommand", Visible 0 }
+];

HelmholtzHABCwithCorners/padeCube.dat

+DefineConstant[
+  dimL = {2.2, Min 1, Step 0.1, Max 5, Name "Input/1Geometry/1Domain length"}
+];
+
+PNT_1_1_1 = 001;
+PNT_2_1_1 = 002;
+PNT_2_1_2 = 003;
+PNT_1_1_2 = 004;
+PNT_1_2_2 = 005;
+PNT_2_2_2 = 006;
+PNT_2_2_1 = 007;
+PNT_1_2_1 = 008;
+
+LIN_0_1_1 = 101;
+LIN_0_1_2 = 102;
+LIN_0_2_2 = 103;
+LIN_0_2_1 = 104;
+LIN_1_0_1 = 105;
+LIN_2_0_1 = 106;
+LIN_2_0_2 = 107;
+LIN_1_0_2 = 108;
+LIN_1_1_0 = 109;
+LIN_1_2_0 = 110;
+LIN_2_2_0 = 111;
+LIN_2_1_0 = 112;
+
+SUR_1_0_0 = 201;
+SUR_2_0_0 = 202;
+SUR_0_1_0 = 203;
+SUR_0_2_0 = 204;
+SUR_0_0_1 = 205;
+SUR_0_0_2 = 206;
+
+SUR_Scatt = 207;
+
+VOL = 301;

HelmholtzHABCwithCorners/padeCube.geo

+Include "padeCube.dat";
+
+/// SPHERE
+
+Point(100) = {     0,     0,     0};
+Point(101) = {-R_SCA,     0,     0};
+Point(102) = { R_SCA,     0,     0};
+Point(103) = {     0,-R_SCA,     0};
+Point(104) = {     0, R_SCA,     0};
+Point(105) = {     0,     0,-R_SCA};
+Point(106) = {     0,     0, R_SCA};
+
+Circle(21) = {102,100,104};
+Circle(22) = {104,100,101};
+Circle(23) = {101,100,103};
+Circle(24) = {103,100,102};
+Circle(25) = {104,100,105};
+Circle(26) = {105,100,103};
+Circle(27) = {103,100,106};
+Circle(28) = {106,100,104};
+Circle(29) = {102,100,106};
+Circle(30) = {106,100,101};
+Circle(31) = {101,100,105};
+Circle(32) = {105,100,102};
+
+Line Loop(41) = { 22, 28,-30};
+Line Loop(42) = { 30, 23, 27};
+Line Loop(43) = {-28,-29, 21};
+Line Loop(44) = {-31,-22, 25};
+Line Loop(45) = {-25,-32,-21};
+Line Loop(46) = {-23, 31, 26};
+Line Loop(47) = {-27, 24, 29};
+Line Loop(48) = {-24, 32,-26};
+Surface(41) = {41};
+Surface(42) = {42};
+Surface(43) = {43};
+Surface(44) = {44};
+Surface(45) = {45};
+Surface(46) = {46};
+Surface(47) = {47};
+Surface(48) = {48};
+
+/// CUBE
+
+X_SCA = dimL/2;
+Y_SCA = dimL/2;
+Z_SCA = dimL/2;
+Point(1) = {     -X_SCA,     -Y_SCA,     -Z_SCA};
+Point(2) = { dimL-X_SCA,     -Y_SCA,     -Z_SCA};
+Point(4) = {     -X_SCA,     -Y_SCA, dimL-Z_SCA};
+Point(3) = { dimL-X_SCA,     -Y_SCA, dimL-Z_SCA};
+Point(5) = {     -X_SCA, dimL-Y_SCA, dimL-Z_SCA};
+Point(6) = { dimL-X_SCA, dimL-Y_SCA, dimL-Z_SCA};
+Point(7) = {     -X_SCA, dimL-Y_SCA,     -Z_SCA};
+Point(8) = { dimL-X_SCA, dimL-Y_SCA,     -Z_SCA};
+
+Line(1) = {1, 2};
+Line(2) = {4, 3};
+Line(3) = {5, 6};
+Line(4) = {7, 8};
+Line(5) = {1, 7};
+Line(6) = {2, 8};
+Line(7) = {3, 6};
+Line(8) = {4, 5};
+Line(9) = {1, 4};
+Line(10) = {7, 5};
+Line(11) = {8, 6};
+Line(12) = {2, 3};
+
+Line Loop(1) = {9, 8, -10, -5};
+Line Loop(2) = {6, 11, -7, -12};
+Line Loop(3) = {1, 12, -2, -9};
+Line Loop(4) = {10, 3, -11, -4};
+Line Loop(5) = {5, 4, -6, -1};
+Line Loop(6) = {2, 7, -3, -8};
+Plane Surface(1) = {1};
+Plane Surface(2) = {2};
+Plane Surface(3) = {3};
+Plane Surface(4) = {4};
+Plane Surface(5) = {5};
+Plane Surface(6) = {6};
+
+/// VOLUME
+
+Surface Loop(1) = {1, 2, 3, 4, 5, 6, -41, -42, -43, -44, -45, -46, -47, -48};
+Volume(1) = {1};
+
+/// PHYSICAL TAGS
+
+Physical Surface(SUR_Scatt) = {-41, -42, -43, -44, -45, -46, -47, -48};
+
+Physical Point(PNT_1_1_1) = {1};
+Physical Point(PNT_2_1_1) = {2};
+Physical Point(PNT_2_1_2) = {3};
+Physical Point(PNT_1_1_2) = {4};
+Physical Point(PNT_1_2_2) = {5};
+Physical Point(PNT_2_2_2) = {6};
+Physical Point(PNT_1_2_1) = {7};
+Physical Point(PNT_2_2_1) = {8};
+
+Physical Line(LIN_0_1_1) = {1};
+Physical Line(LIN_0_1_2) = {2};
+Physical Line(LIN_0_2_2) = {3};
+Physical Line(LIN_0_2_1) = {4};
+Physical Line(LIN_1_0_1) = {5};
+Physical Line(LIN_2_0_1) = {6};
+Physical Line(LIN_2_0_2) = {7};
+Physical Line(LIN_1_0_2) = {8};
+Physical Line(LIN_1_1_0) = {9};
+Physical Line(LIN_1_2_0) = {10};
+Physical Line(LIN_2_2_0) = {11};
+Physical Line(LIN_2_1_0) = {12};
+
+Physical Surface(SUR_1_0_0) = {1};
+Physical Surface(SUR_2_0_0) = {2};
+Physical Surface(SUR_0_1_0) = {3};
+Physical Surface(SUR_0_2_0) = {4};
+Physical Surface(SUR_0_0_1) = {5};
+Physical Surface(SUR_0_0_2) = {6};
+
+Physical Volume(VOL) = {1};

HelmholtzHABCwithCorners/padeDisk.dat

+BND_Neumann        = 0;
+BND_Sommerfeld     = 1;
+BND_Second         = 2;
+BND_Pade           = 3;
+BND_CRBC           = 4;
+BND_PML            = 5;
+
+DefineConstant[
+  BND_TYPE = {BND_PML,
+    Name "Input/5Model/02Boundary condition (edges)", Highlight "Blue",
+    Choices {BND_Neumann    = "Homogeneous Neumann",
+             BND_Sommerfeld = "Sommerfeld ABC",
+             BND_Second     = "Second-order ABC",
+             BND_Pade       = "Pade ABC",
+             BND_CRBC       = "CRBC",
+             BND_PML        = "PML"}}
+];
+
+DefineConstant[
+  R_DOM  = { 1.1, Min 1, Step 0.1, Max 10,   Name "Input/1Geometry/1Domain length"},
+  Npml   = { 5,   Min 1, Step 1,   Max 5,    Name "Input/5Model/03PML: Thickness (N*Lc)", Visible (BND_TYPE == BND_PML)},
+  rotPml = { 15,  Min 0, Step 2.5, Max 92.5, Name "Input/5Model/03PML: Rotation", Visible (BND_TYPE == BND_PML)}
+];
+
+Lpml = LC*Npml;
+
+BND_Scatt = 201;
+BND_Dom = 202;
+BND_PmlExt = 203;
+
+DOM = 301;
+DOM_PML = 302;

HelmholtzHABCwithCorners/padeDisk.geo

+Include "padeDisk.dat";
+
+Point(0) = {0, 0, 0};
+
+Point(1) = {-R_SCA,     0, 0};
+Point(2) = {     0,-R_SCA, 0};
+Point(3) = { R_SCA,     0, 0};
+Point(4) = {     0, R_SCA, 0};
+
+Circle(1) = {1, 0, 2};
+Circle(2) = {2, 0, 3};
+Circle(3) = {3, 0, 4};
+Circle(4) = {4, 0, 1};
+
+If((R_DOM-R_SCA) > LC)
+  Point(5) = {-R_DOM,     0, 0};
+  Point(6) = {     0,-R_DOM, 0};
+  Point(7) = { R_DOM,     0, 0};
+  Point(8) = {     0, R_DOM, 0};
+
+  Circle(5) = {5, 0, 6};
+  Circle(6) = {6, 0, 7};
+  Circle(7) = {7, 0, 8};
+  Circle(8) = {8, 0, 5};
+
+  Line Loop(1) = {1, 2, 3, 4, -5, -6, -7, -8};
+  Plane Surface(1) = {1};
+
+  If(BND_TYPE == BND_PML)
+  If(Npml > 0)
+    Extrude { Line{5, 6, 7, 8}; Layers{Npml,Npml*LC}; Recombine; }
+    Physical Line(BND_PmlExt) = {9, 13, 17, 21};
+    Physical Surface(DOM_PML) = {12, 16, 20, 24};
+  Else
+    Physical Line(BND_PmlExt) = {5, 6, 7, 8};
+    Physical Surface(DOM_PML) = {};
+  EndIf
+  EndIf
+
+  Physical Line(BND_Dom) = {5, 6, 7, 8};
+  Physical Surface(DOM) = {1};
+Else
+
+  If(BND_TYPE == BND_PML)
+    Extrude { Line{1, 2, 3, 4}; Layers{Npml,Npml*LC}; Recombine; }
+    Physical Line(BND_PmlExt) = {5, 9, 13, 17};
+    Physical Surface(DOM_PML) = {8, 12, 16, 20};
+  EndIf
+
+  Physical Line(BND_Dom) = {};
+  Physical Surface(DOM) = {};
+EndIf
+
+Physical Line(BND_Scatt) = {1, 2, 3, 4};

HelmholtzHABCwithCorners/padeDisk.pro

+Include "padeDisk.dat";
+
+//==================================================================================================
+// OPTIONS and PARAMETERS
+//==================================================================================================
+
+DefineConstant[
+  nPade = {4, Choices {0, 1, 2, 3, 4, 5, 6},
+    Name "Input/5Model/05Pade: Number of fields",
+    Visible (BND_TYPE == BND_Pade)},
+  thetaPadeInput = {3, Choices {0, 1, 2, 3, 4},
+    Name "Input/5Model/06Pade: Rotation of branch cut",
+    Visible (BND_TYPE == BND_Pade)}
+];
+
+If(BND_TYPE == BND_Pade)
+  If(thetaPadeInput == 0)
+    thetaPade = 0;
+  ElseIf(thetaPadeInput == 1)
+    thetaPade = Pi/8;
+  ElseIf(thetaPadeInput == 2)
+    thetaPade = Pi/4;
+  ElseIf(thetaPadeInput == 3)
+    thetaPade = Pi/3;
+  ElseIf(thetaPadeInput == 4)
+    thetaPade = Pi/2;
+  EndIf
+  mPade = 2*nPade+1;
+  For j In{1:nPade}
+    cPade~{j} = Tan[j*Pi/mPade]^2;
+  EndFor
+Else
+  nPade = 0;
+  thetaPadeInput = 0;
+EndIf
+If(BND_TYPE == BND_PML)
+  nPade = Npml;
+EndIf
+
+Group {
+  Dom    = Region[{DOM}];
+  BndSca = Region[{BND_Scatt}];
+If(BND_TYPE == BND_PML)
+  DomPml = Region[{DOM_PML}];
+  BndExt = Region[{BND_PmlExt}];
+Else
+  DomPml = Region[{}];
+  BndExt = Region[{BND_Dom}];
+EndIf
+  DomAll = Region[{Dom,DomPml,BndSca,BndExt}];
+}
+
+Function {
+
+If(BND_TYPE == BND_Pade)
+  kEps[] = WAVENUMBER + I[] * 0.39 * WAVENUMBER^(1/3) * (1/R_DOM)^(2/3);
+EndIf
+
+If(BND_TYPE == BND_Pade)
+  ExpPTheta[]  = Complex[Cos[ thetaPade],Sin[ thetaPade]];
+  ExpMTheta[]  = Complex[Cos[-thetaPade],Sin[-thetaPade]];
+  ExpPTheta2[] = Complex[Cos[thetaPade/2.],Sin[thetaPade/2.]];
+  For i In{1:nPade}
+  For j In{1:nPade}
+    coefA~{i}~{j}[] = 2./mPade * cPade~{j} * (cPade~{i}-1+ExpMTheta[]) / (cPade~{i}+cPade~{j}+ExpMTheta[]);
+    coefB~{i}~{j}[] = 2./mPade * cPade~{j} * (-1-cPade~{i}) / (cPade~{i}+cPade~{j}+ExpMTheta[]);
+  EndFor
+  EndFor
+EndIf
+
+If(BND_TYPE == BND_PML)
+  rLoc[DomPml] = R[]-R_DOM;
+  absFuncS[DomPml] = 1/(Lpml-rLoc[]);
+  absFuncF[DomPml] = -Log[1-rLoc[]/Lpml];
+  //absFuncS[DomPml] = 1/(Lpml-rLoc[]) - 1/Lpml;
+  //absFuncF[DomPml] = -Log[1-rLoc[]/Lpml] - rLoc[]/Lpml;
+  If(rotPml < 91)
+    rot[DomPml] = Complex[Sin[rotPml*Pi/180.], Cos[rotPml*Pi/180.]]; // I (rotPml=0, prop) - 1 (rotPml=Pi/2, evan)
+  Else
+    rot[DomPml] = Complex[1., 1.];
+  EndIf
+  s1[DomPml] = 1 + rot[] * absFuncS[]/k[];
+  s2[DomPml] = 1 + rot[] * (1/R[]) * absFuncF[]/k[];
+  nVec[DomPml] = XYZ[]/R[];
+  tVec[DomPml] = nVec[] /\ Vector[0,0,1];
+  nTen[DomPml] = Tensor[ CompX[nVec[]]*CompX[nVec[]], CompX[nVec[]]*CompY[nVec[]], CompX[nVec[]]*CompZ[nVec[]],
+                         CompY[nVec[]]*CompX[nVec[]], CompY[nVec[]]*CompY[nVec[]], CompY[nVec[]]*CompZ[nVec[]],
+                         CompZ[nVec[]]*CompX[nVec[]], CompZ[nVec[]]*CompY[nVec[]], CompZ[nVec[]]*CompZ[nVec[]] ];
+  tTen[DomPml] = Tensor[ CompX[tVec[]]*CompX[tVec[]], CompX[tVec[]]*CompY[tVec[]], CompX[tVec[]]*CompZ[tVec[]],
+                         CompY[tVec[]]*CompX[tVec[]], CompY[tVec[]]*CompY[tVec[]], CompY[tVec[]]*CompZ[tVec[]],
+                         CompZ[tVec[]]*CompX[tVec[]], CompZ[tVec[]]*CompY[tVec[]], CompZ[tVec[]]*CompZ[tVec[]] ];
+  pmlScal[DomPml] = s1[]*s2[];
+  pmlTens[DomPml] = (s2[]/s1[]) * nTen[] + (s1[]/s2[]) * tTen[];
+EndIf
+}
+
+//==================================================================================================
+// FONCTION SPACES with CONSTRAINTS
+//==================================================================================================
+
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+Constraint {
+  { Name DirichletBC; Case {
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+    { Region BndSca; Value -f_inc[]; }
+EndIf
+  }}
+}
+EndIf
+
+FunctionSpace {
+  { Name H_ref; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{DomAll}]; Entity NodesOf[All]; }}
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+    Constraint {{ NameOfCoef pi; EntityType NodesOf; NameOfConstraint DirichletBC; }}
+EndIf
+  }
+  { Name H_num; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{DomAll}]; Entity NodesOf[All]; }}
+If((FLAG_SIGNAL_BC == SIGNAL_Dirichlet) || (BND_TYPE == BND_PML))
+    Constraint {{ NameOfCoef pi; EntityType NodesOf; NameOfConstraint DirichletBC; }}
+EndIf
+  }
+If(BND_TYPE == BND_Pade)
+  For i In {1:nPade}
+  { Name H~{i}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
+  EndFor
+EndIf
+}
+
+//==================================================================================================
+// FORMULATIONS
+//==================================================================================================
+
+Formulation {
+
+  { Name NumSol; Type FemEquation;
+    Quantity {
+      { Name u_num; Type Local; NameOfSpace H_num; }
+If(BND_TYPE == BND_Pade)
+    For i In {1:nPade}
+      { Name u~{i}; Type Local; NameOfSpace H~{i}; }
+    EndFor
+EndIf
+    }
+    Equation {
+
+// Helmholtz
+
+      Galerkin{ [ Dof{d u_num}    , {d u_num} ]; In Dom;    Jacobian JVol; Integration I1; }
+      Galerkin{ [ -k[]^2*Dof{u_num} , {u_num} ]; In Dom;    Jacobian JVol; Integration I1; }
+If(FLAG_SIGNAL_BC == SIGNAL_Neumann)
+      Galerkin{ [ -df_inc[]         , {u_num} ]; In BndSca; Jacobian JSur; Integration I1; }
+EndIf
+If(BND_TYPE == BND_PML)
+      Galerkin{ [ pmlTens[] * Dof{d u_num}      , {d u_num} ]; In DomPml; Jacobian JVol; Integration I1; }
+      Galerkin{ [ -k[]^2 * pmlScal[] * Dof{u_num} , {u_num} ]; In DomPml; Jacobian JVol; Integration I1; }
+EndIf
+
+// Sommerfeld ABC
+
+If(BND_TYPE == BND_Sommerfeld)
+      Galerkin{ [ -I[]*k[]*Dof{u_num} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+EndIf
+
+// Second-order ABC
+
+If(BND_TYPE == BND_Second)
+      Galerkin { [ - I[]*k[] * Dof{u_num}           , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+      Galerkin { [ - 1/(2*I[]*k[]) * Dof{d u_num} , {d u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+EndIf
+
+// HABC
+
+If(BND_TYPE == BND_Pade)
+      Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+    For i In{1:nPade}
+      Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u~{i}} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+      Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u_num} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+
+      Galerkin { [ Dof{d u~{i}}                                   , {d u~{i}} ]; In BndExt; Jacobian JSur; Integration I1; }
+      Galerkin { [ -kEps[]^2*(ExpPTheta[]*cPade~{i}+1) * Dof{u~{i}} , {u~{i}} ]; In BndExt; Jacobian JSur; Integration I1; }
+      Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{i}+1) * Dof{u_num} , {u~{i}} ]; In BndExt; Jacobian JSur; Integration I1; }
+    EndFor
+EndIf
+
+    }
+  }
+
+  { Name ProjSol; Type FemEquation;
+    Quantity {
+      { Name u_refProj; Type Local; NameOfSpace H_ref; }
+    }
+    Equation {
+      Galerkin{ [ Dof{u_refProj} , {u_refProj} ]; In Dom; Jacobian JVol; Integration I1; }
+      Galerkin{ [ -f_ref[]       , {u_refProj} ]; In Dom; Jacobian JVol; Integration I1; }
+    }
+  }
+}
+
+//==================================================================================================
+// RESOLUTION
+//==================================================================================================
+
+Resolution {
+  { Name NumSol;
+    System {
+      { Name A; NameOfFormulation NumSol; Type Complex; }
+    }
+    Operation {
+      Generate[A]; Solve[A]; SaveSolution[A];
+    }
+  }
+  { Name ProjSol;
+    System {
+      { Name A; NameOfFormulation ProjSol; Type Complex; }
+    }
+    Operation {
+      Generate[A]; Solve[A]; SaveSolution[A];
+    }
+  }
+}
+
+//==================================================================================================
+// POSTPRO / POSTOP
+//==================================================================================================
+
+PostProcessing {
+  { Name NumSol; NameOfFormulation NumSol;
+    Quantity {
+//      { Name param1; Value { Local { [ rLoc[] ]; In DomPml; Jacobian JVol; }}}
+//      { Name param2; Value { Local { [ nVec[] ]; In DomPml; Jacobian JVol; }}}
+//      { Name param3; Value { Local { [ tVec[] ]; In DomPml; Jacobian JVol; }}}
+      { Name u_ref; Value { Local { [ f_ref[] ]; In Region[{Dom,BndSca}]; Jacobian JVol; }}}
+      { Name u_num~{BND_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ {u_num} ];         In Region[{Dom,DomPml}]; Jacobian JVol; }}}
+      { Name u_err~{BND_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ f_ref[]-{u_num} ]; In Region[{Dom,BndSca}]; Jacobian JVol; }}}
+    }
+  }
+  { Name Errors; NameOfFormulation NumSol;
+    Quantity {
+      { Name energy;       Value { Integral { [ Abs[f_ref[]]^2 ];                          In Dom; Jacobian JVol; Integration I1; }}}
+      { Name error2;       Value { Integral { [ Abs[f_ref[]-{u_num}]^2 ];                  In Dom; Jacobian JVol; Integration I1; }}}
+      { Name errorAbs;     Value { Term { Type Global; [Sqrt[$error2Var]] ;                In Dom; Jacobian JVol; }}}
+      { Name errorRel;     Value { Term { Type Global; [Sqrt[$error2Var/$energyVar]] ;     In Dom; Jacobian JVol; }}}
+    }
+  }
+  { Name ProjError; NameOfFormulation ProjSol;
+    Quantity {
+      { Name energy;       Value { Integral { [ Abs[f_ref[]]^2 ];                          In Dom; Jacobian JVol; Integration I1; }}}
+      { Name errorProj2;   Value { Integral { [ Abs[f_ref[]-{u_refProj}]^2 ];              In Dom; Jacobian JVol; Integration I1; }}}
+      { Name errorProjAbs; Value { Term { Type Global; [Sqrt[$errorProj2Var]] ;            In Dom; Jacobian JVol; }}}
+      { Name errorProjRel; Value { Term { Type Global; [Sqrt[$errorProj2Var/$energyVar]] ; In Dom; Jacobian JVol; }}}
+    }
+  }
+  { Name DtNError; NameOfFormulation NumSol;
+    Quantity {
+      { Name energyDtN;       Value { Integral { [ Abs[f_ref[]]^2 ];                            In BndSca; Jacobian JLin; Integration I1; }}}
+      { Name errorDtN2;       Value { Integral { [ Abs[f_ref[]-{u_num}]^2 ];                    In BndSca; Jacobian JLin; Integration I1; }}}
+      { Name normDtN;         Value { Term { Type Global; [Sqrt[$energyDtNVar]] ;               In BndSca; Jacobian JLin; }}}
+      { Name errorDtNAbs;     Value { Term { Type Global; [Sqrt[$errorDtN2Var]] ;               In BndSca; Jacobian JLin; }}}
+      { Name errorDtNRel;     Value { Term { Type Global; [Sqrt[$errorDtN2Var/$energyDtNVar]] ; In BndSca; Jacobian JLin; }}}
+    }
+  }
+}
+
+PostOperation{
+  { Name NumSol; NameOfPostProcessing NumSol;
+    Operation {
+//      Print [param1, OnElementsOf DomPml, File "out/param1.pos"];
+//      Print [param2, OnElementsOf DomPml, File "out/param2.pos"];
+//      Print [param3, OnElementsOf DomPml, File "out/param3.pos"];
+      tmp1 = Sprintf("out/solNum_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
+      tmp2 = Sprintf("out/solErr_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
+      Print [u_ref, OnElementsOf Region[{DOM,BND_Scatt}], File "out/solRef.pos"];
+      Print [u_num~{BND_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Region[{Dom,DomPml}], File tmp1];
+      Print [u_err~{BND_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Region[{Dom,BndSca}], File tmp2];
+    }
+  }
+  { Name Errors; NameOfPostProcessing Errors;
+    Operation {
+      tmp3 = Sprintf("out/errorAbs_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
+      tmp4 = Sprintf("out/errorRel_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
+      Print [energy[Dom], OnRegion Dom, Format Table, StoreInVariable $energyVar];
+      Print [error2[Dom], OnRegion Dom, Format Table, StoreInVariable $error2Var];
+      Print [errorAbs,    OnRegion Dom, Format Table, SendToServer "Output/1L2-Error (absolute)", File > tmp3];
+      Print [errorRel,    OnRegion Dom, Format Table, SendToServer "Output/2L2-Error (relative)", File > tmp4];
+    }
+  }
+  { Name ProjError; NameOfPostProcessing ProjError;
+    Operation {
+      tmp5 = Sprintf("out/errorProjAbs_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
+      tmp6 = Sprintf("out/errorProjRel_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
+      Print [energy[Dom],     OnRegion Dom, Format Table, StoreInVariable $energyVar];
+      Print [errorProj2[Dom], OnRegion Dom, Format Table, StoreInVariable $errorProj2Var];
+      Print [errorProjAbs,    OnRegion Dom, Format Table, SendToServer "Output/3L2-ErrorProj (absolute)", File > tmp5];
+      Print [errorProjRel,    OnRegion Dom, Format Table, SendToServer "Output/4L2-ErrorProj (relative)", File > tmp6];
+    }
+  }
+  { Name DtNError; NameOfPostProcessing DtNError;
+    Operation {
+      tmp7 = Sprintf("out/normDtNAbs_%g_%g_%g_%g.dat",  FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
+      tmp8 = Sprintf("out/errorDtNAbs_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
+      tmp9 = Sprintf("out/errorDtNRel_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
+      Print [energyDtN[BndSca], OnRegion BndSca, Format Table, StoreInVariable $energyDtNVar];
+      Print [errorDtN2[BndSca], OnRegion BndSca, Format Table, StoreInVariable $errorDtN2Var];
+      Print [normDtN,           OnRegion BndSca, Format Table, SendToServer "Output/4L2-DtN-Norm", File > tmp7];
+      Print [errorDtNAbs,       OnRegion BndSca, Format Table, SendToServer "Output/5L2-DtN-Error (absolute)", File > tmp8];
+      Print [errorDtNRel,       OnRegion BndSca, Format Table, SendToServer "Output/6L2-DtN-Error (relative)", File > tmp9];
+    }
+  }
+}
+
+DefineConstant[
+  R_ = {"NumSol", Name "GetDP/1ResolutionChoices",    Visible 1, Choices {"NumSol", "ProjSol"}},
+  P_ = {"NumSol", Name "GetDP/2PostOperationChoices", Visible 1, Choices {"NumSol", "Errors", "ProjError", "DtNError"}}
+];

HelmholtzHABCwithCorners/padePieWedge.dat

+REF_Ana = 0;
+REF_Num = 1;
+
+DefineConstant[
+  dimL     = { 3.3, Min 0.1, Step 0.1, Max  20, Name "Input/1Geometry/1Wedge length"},
+  alphaDom = {  90, Min 60,  Step 10 , Max 180, Name "Input/1Geometry/2Wedge angle"},
+  distSca  = { 2.2, Min 0.1, Step 0.1, Max  10, Name "Input/1Geometry/3Scatterer distance"},
+//  X_SCA = { 1.1, Min 1, Step 0.01, Max 10, Name "Input/1Geometry/2Scatterer position (x0)"},
+//  Y_SCA = { 1.1, Min 1, Step 0.01, Max 10, Name "Input/1Geometry/3Scatterer position (y0)"}
+  FLAG_REF = {REF_Num,
+    Name "Input/1Reference solution",
+    Choices {REF_Ana = "Analytic solution", REF_Num = "Numeric solution"}}
+];
+
+alphaDomSave = alphaDom;
+alphaDom = alphaDom*Pi/180;
+alphaSca = alphaDom/2;
+//dimL = (R_SCA*1.1)*(1.+1./Sin[alphaSca]);
+//distSca = (R_SCA*1.1)/Sin[alphaSca];
+X_SCA = distSca * Sin[alphaSca];
+Y_SCA = distSca * Cos[alphaSca];
+
+X~{1}   =                      - X_SCA;
+Y~{1}   =                 dimL - Y_SCA;
+X~{100} = Sin[alphaDom/2]*dimL - X_SCA;
+Y~{100} = Cos[alphaDom/2]*dimL - Y_SCA;
+X~{2}   = Sin[alphaDom  ]*dimL - X_SCA;
+Y~{2}   = Cos[alphaDom  ]*dimL - Y_SCA;
+X~{3}   =                      - X_SCA;
+Y~{3}   =                      - Y_SCA;
+X~{200} =-Sin[alphaDom/2]*dimL - X_SCA;
+Y~{200} =-Cos[alphaDom/2]*dimL - Y_SCA;
+
+CRN_1 = 101;
+CRN_2 = 102;
+CRN_3 = 103;
+
+BND_1 = 201;
+BND_2 = 202;
+BND_3 = 203;
+
+BND_Scatt = 205;
+BND_EXT   = 206;
+
+DOM = 301;
+DOM_EXT = 302;

HelmholtzHABCwithCorners/padePieWedge.geo

+Include "padePieWedge.dat";
+
+Point(0) = {0, 0, 0};
+
+Point(1) = {X~{1}, Y~{1}, 0};
+Point(2) = {X~{2}, Y~{2}, 0};
+Point(3) = {X~{3}, Y~{3}, 0};
+Point(100) = {X~{100}, Y~{100}, 0};
+
+Point(5) = {-R_SCA,     0, 0};
+Point(6) = {     0,-R_SCA, 0};
+Point(7) = { R_SCA,     0, 0};
+Point(8) = {     0, R_SCA, 0};
+
+Circle(0) = {2, 3, 100};
+Circle(1) = {100, 3, 1};
+Line(2) = {3, 2};
+Line(3) = {1, 3};
+
+Circle(5) = {5, 0, 6};
+Circle(6) = {6, 0, 7};
+Circle(7) = {7, 0, 8};
+Circle(8) = {8, 0, 5};
+
+Line Loop(1) = {0, 1, 2, 3, -5, -6, -7, -8};
+Plane Surface(1) = {1};
+
+If(FLAG_REF == REF_Num)
+  Point(200) = {X~{200}, Y~{200}, 0};
+  Circle(10) = {1, 3, 200};
+  Circle(11) = {200, 3, 2};
+  Line Loop(2) = {-2, -3, 10, 11};
+  Plane Surface(2) = {2};
+EndIf
+
+Physical Point(CRN_1) = {1}; // Hybrid1
+Physical Point(CRN_2) = {2}; // Hybrid2
+Physical Point(CRN_3) = {3}; // Pade
+
+Physical Line(BND_1) = {0, 1}; // BGT
+Physical Line(BND_2) = {2};    // Pade-Left
+Physical Line(BND_3) = {3};    // Pade-Down
+Physical Line(BND_Scatt) = {5,6,7,8};
+
+SetOrder 1;
+Mesh.ElementOrder = 1;
+Mesh 1;
+Save "mainCurv.msh";
+SetOrder ORDER;
+Mesh.ElementOrder = ORDER;
+
+Physical Surface(DOM) = {1};
+
+If(FLAG_REF == REF_Num)
+  Physical Line(BND_EXT) = {10, 11};
+  Physical Surface(DOM_EXT) = {2};
+EndIf
+
+

HelmholtzHABCwithCorners/padePolygon.dat

+REF_Ana = 1;
+REF_Num = 2;
+
+DefineConstant[
+  DOM_midR       = { 1.65, Min 1, Step 0.1, Max 10,         Name "Input/1Geometry/1Polygon: mid-radius"},
+  DOM_Num        = { 3, Choices {3, 4, 5, 6, 8, 12, 18, 1}, Name "Input/1Geometry/2Polygon: number edges"},
+  DOM_AngleEdges = { 180.-360./DOM_Num,                     Name "Input/1Geometry/32Polygon: angle edges", ReadOnly 1},
+  FLAG_REF = {REF_Ana,
+    Name "Input/1Reference solution",
+    Choices {REF_Ana = "Analytic solution", REF_Num = "Numeric solution"}},
+  L_EXT  = { 5, Min 1, Step 0.1, Max 10, Name "Input/2Geometry/0Reference length",
+    Visible (FLAG_REF == REF_Num)}
+];
+
+angleInterior = 360/DOM_Num;
+angleEdges = 180-angleInterior;
+
+DOM_R = 2./(1.+Cos[Pi/DOM_Num])*DOM_midR;
+If(DOM_Num == 1)
+  DOM_R = DOM_midR;
+EndIf
+
+For i In{1:DOM_Num}
+  X~{i} = DOM_R*Cos[(i-0.5)*angleInterior*Pi/180-Pi/2];
+  Y~{i} = DOM_R*Sin[(i-0.5)*angleInterior*Pi/180-Pi/2];
+  CRN~{i} = 100+i;
+  BND~{i} = 200+i;
+EndFor
+BND_SCA = 299;
+DOM     = 301;
+
+For i In{1:4}
+  CRN_EXT~{i} = 400+i;
+  BND_EXT~{i} = 500+i;
+EndFor
+BND_EXT = 599;
+DOM_EXT = 601;

HelmholtzHABCwithCorners/padePolygon.geo

+Include "padePolygon.dat";
+
+Point(0) = {0, 0, 0};
+Point(1) = {-R_SCA,     0, 0};
+Point(2) = {     0,-R_SCA, 0};
+Point(3) = { R_SCA,     0, 0};
+Point(4) = {     0, R_SCA, 0};
+Circle(1) = {1, 0, 2};
+Circle(2) = {2, 0, 3};
+Circle(3) = {3, 0, 4};
+Circle(4) = {4, 0, 1};
+llScat[] = {-1, -2, -3, -4};
+
+If(DOM_Num == 1)
+  Point(11) = {-DOM_R,     0, 0};
+  Point(12) = {     0,-DOM_R, 0};
+  Point(13) = { DOM_R,     0, 0};
+  Point(14) = {     0, DOM_R, 0};
+  Circle(11) = {11, 0, 12};
+  Circle(12) = {12, 0, 13};
+  Circle(13) = {13, 0, 14};
+  Circle(14) = {14, 0, 11};
+  llBnd[] = {-11, -12, -13, -14};
+Else
+For iCrn In{1:DOM_Num}
+  pBnd~{iCrn} = newp;
+  Point(pBnd~{iCrn}) = {X~{iCrn}, Y~{iCrn}, 0};
+EndFor
+For iEdg In{1:DOM_Num}
+  iCrn1 = (iEdg == 1) ? DOM_Num : iEdg-1;
+  iCrn2 = iEdg;
+  lBnd~{iEdg} = newl;
+  Line(lBnd~{iEdg}) = {pBnd~{iCrn1}, pBnd~{iCrn2}};
+  llBnd[] += lBnd~{iEdg};
+EndFor
+EndIf
+
+Line Loop(1) = {llBnd[], llScat[]};
+Plane Surface(1) = {1};
+
+If(DOM_Num == 1)
+  Physical Line(BND~{1}) = {11, 12, 13, 14};
+Else
+For i In{1:DOM_Num}
+  Physical Point(CRN~{i}) = {pBnd~{i}};
+  Physical Line(BND~{i})  = {lBnd~{i}};
+EndFor
+EndIf
+
+SetOrder 1;
+Mesh.ElementOrder = 1;
+Mesh 1;
+Save "mainCurv.msh";
+SetOrder ORDER;
+Mesh.ElementOrder = ORDER;
+
+Physical Line(BND_SCA) = {1,2,3,4};
+Physical Surface(DOM) = {1};
+
+If(FLAG_REF == REF_Num)
+  pExt~{1} = newp; Point(pExt~{1}) = {-L_EXT/2,-L_EXT/2, 0};
+  pExt~{2} = newp; Point(pExt~{2}) = { L_EXT/2,-L_EXT/2, 0};
+  pExt~{3} = newp; Point(pExt~{3}) = { L_EXT/2, L_EXT/2, 0};
+  pExt~{4} = newp; Point(pExt~{4}) = {-L_EXT/2, L_EXT/2, 0};
+  lExt~{1} = newl; Line(lExt~{1}) = {pExt~{4}, pExt~{1}};
+  lExt~{2} = newl; Line(lExt~{2}) = {pExt~{1}, pExt~{2}};
+  lExt~{3} = newl; Line(lExt~{3}) = {pExt~{2}, pExt~{3}};
+  lExt~{4} = newl; Line(lExt~{4}) = {pExt~{3}, pExt~{4}};
+  For i In{1:4}
+    Physical Point(CRN_EXT~{i}) = {pExt~{i}};
+    Physical Line(BND_EXT~{i}) = {lExt~{i}};
+    llExt[] += lExt~{i};
+  EndFor
+  Line Loop(2) = {llBnd[], llExt[]};
+  Plane Surface(2) = {2};
+  Physical Surface(DOM_EXT) = {2};
+EndIf

HelmholtzHABCwithCorners/padePolyhedron.dat

+DefineConstant[
+  dimRmean  = {2, Min 1, Step 0.1, Max 5, Name "Input/1Geometry/3Polyhedron: midradius"},
+  FAC_Num = {6,
+    Name "Input/1Geometry/2Polyhedron: number faces", Highlight "Blue",
+    Choices {4  = "Tetrahedron",
+             6  = "Cube",
+             8  = "Octahedron",
+             12 = "Dodecahedron",
+             20 = "Icosahedron",
+             1  = "Sphere"} }
+];
+
+GoldRatio = 1.61803398875;
+If(FAC_Num == 4)
+  dimRext = dimRmean * Sin[Pi/4]/Cos[Pi/3] * Tan[Pi/3]/Sqrt[2];
+  dimRint = dimRext/(Tan[Pi/3]*Tan[Pi/3]);
+ElseIf(FAC_Num == 6)
+  dimRext = dimRmean * Sin[Pi/6]/Cos[Pi/4] * Tan[Pi/3];
+  dimRint = dimRext/(Tan[Pi/4]*Tan[Pi/3]);
+ElseIf(FAC_Num == 8)
+  dimRext = dimRmean * Sin[Pi/6]/Cos[Pi/3] * Tan[Pi/4]*Sqrt[2];
+  dimRint = dimRext/(Tan[Pi/3]*Tan[Pi/4]);
+ElseIf(FAC_Num == 12)
+  dimRext = dimRmean * Sin[Pi/10]/Cos[Pi/5] * Tan[Pi/3] * GoldRatio;
+  dimRint = dimRext/(Tan[Pi/5]*Tan[Pi/3]);
+ElseIf(FAC_Num == 20)
+  dimRext = dimRmean * Sin[Pi/10]/Cos[Pi/3] * Tan[Pi/5] * GoldRatio*GoldRatio;
+  dimRint = dimRext/(Tan[Pi/3]*Tan[Pi/5]);
+ElseIf(FAC_Num == 1)
+  dimRext = dimRmean;
+  dimRint = dimRmean;
+EndIf
+Printf("... Polyhedron circumradius =", dimRext);
+Printf("... Polyhedron inradius     =", dimRint);
+
+// ADJACENCY
+
+If(FAC_Num == 4)
+  EDG_Num = 6;
+  EdgNeigh~{1}~{1} = 1; EdgNeigh~{1}~{2} = 2; // Neighboring faces
+  EdgNeigh~{2}~{1} = 1; EdgNeigh~{2}~{2} = 3;
+  EdgNeigh~{3}~{1} = 2; EdgNeigh~{3}~{2} = 3;
+  EdgNeigh~{4}~{1} = 1; EdgNeigh~{4}~{2} = 4;
+  EdgNeigh~{5}~{1} = 2; EdgNeigh~{5}~{2} = 4;
+  EdgNeigh~{6}~{1} = 3; EdgNeigh~{6}~{2} = 4;
+  CRN_Num = 4;
+  CRN_NumEdgPerCrn = 3;
+  CrnNeigh~{1}~{1} = 1; CrnNeigh~{1}~{2} = 2; CrnNeigh~{1}~{3} = 3; // Neighboring edges
+  CrnNeigh~{2}~{1} = 1; CrnNeigh~{2}~{2} = 4; CrnNeigh~{2}~{3} = 5;
+  CrnNeigh~{3}~{1} = 2; CrnNeigh~{3}~{2} = 4; CrnNeigh~{3}~{3} = 6;
+  CrnNeigh~{4}~{1} = 3; CrnNeigh~{4}~{2} = 5; CrnNeigh~{4}~{3} = 6;
+EndIf
+
+If(FAC_Num == 6)
+  EDG_Num = 12;
+  EdgNeigh~{1}~{1} = 3; EdgNeigh~{1}~{2} = 5;
+  EdgNeigh~{2}~{1} = 3; EdgNeigh~{2}~{2} = 6;
+  EdgNeigh~{3}~{1} = 4; EdgNeigh~{3}~{2} = 6;
+  EdgNeigh~{4}~{1} = 4; EdgNeigh~{4}~{2} = 5;
+  EdgNeigh~{5}~{1} = 1; EdgNeigh~{5}~{2} = 5;
+  EdgNeigh~{6}~{1} = 2; EdgNeigh~{6}~{2} = 5;
+  EdgNeigh~{7}~{1} = 2; EdgNeigh~{7}~{2} = 6;
+  EdgNeigh~{8}~{1} = 1; EdgNeigh~{8}~{2} = 6;
+  EdgNeigh~{9}~{1} = 1; EdgNeigh~{9}~{2} = 3;
+  EdgNeigh~{10}~{1} = 1; EdgNeigh~{10}~{2} = 4;
+  EdgNeigh~{11}~{1} = 2; EdgNeigh~{11}~{2} = 4;
+  EdgNeigh~{12}~{1} = 2; EdgNeigh~{12}~{2} = 3;
+  CRN_Num = 8;
+  CRN_NumEdgPerCrn = 3;
+  CrnNeigh~{1}~{1} =  9; CrnNeigh~{1}~{2} = 5; CrnNeigh~{1}~{3} = 1;
+  CrnNeigh~{2}~{1} = 12; CrnNeigh~{2}~{2} = 6; CrnNeigh~{2}~{3} = 1;
+  CrnNeigh~{3}~{1} = 12; CrnNeigh~{3}~{2} = 7; CrnNeigh~{3}~{3} = 2;
+  CrnNeigh~{4}~{1} =  9; CrnNeigh~{4}~{2} = 8; CrnNeigh~{4}~{3} = 2;
+  CrnNeigh~{5}~{1} = 10; CrnNeigh~{5}~{2} = 8; CrnNeigh~{5}~{3} = 3;
+  CrnNeigh~{6}~{1} = 11; CrnNeigh~{6}~{2} = 7; CrnNeigh~{6}~{3} = 3;
+  CrnNeigh~{7}~{1} = 10; CrnNeigh~{7}~{2} = 5; CrnNeigh~{7}~{3} = 4;
+  CrnNeigh~{8}~{1} = 11; CrnNeigh~{8}~{2} = 6; CrnNeigh~{8}~{3} = 4;
+EndIf
+
+If(FAC_Num == 8)
+  EDG_Num = 12;
+  EdgNeigh~{1}~{1} = 1; EdgNeigh~{1}~{2} = 4;
+  EdgNeigh~{2}~{1} = 1; EdgNeigh~{2}~{2} = 2;
+  EdgNeigh~{3}~{1} = 2; EdgNeigh~{3}~{2} = 3;
+  EdgNeigh~{4}~{1} = 3; EdgNeigh~{4}~{2} = 4;
+  EdgNeigh~{5}~{1} = 1; EdgNeigh~{5}~{2} = 5;
+  EdgNeigh~{6}~{1} = 2; EdgNeigh~{6}~{2} = 6;
+  EdgNeigh~{7}~{1} = 3; EdgNeigh~{7}~{2} = 7;
+  EdgNeigh~{8}~{1} = 4; EdgNeigh~{8}~{2} = 8;
+  EdgNeigh~{9}~{1} = 5; EdgNeigh~{9}~{2} = 8;
+  EdgNeigh~{10}~{1} = 5; EdgNeigh~{10}~{2} = 6;
+  EdgNeigh~{11}~{1} = 6; EdgNeigh~{11}~{2} = 7;
+  EdgNeigh~{12}~{1} = 7; EdgNeigh~{12}~{2} = 8;
+  CRN_Num = 6;
+  CRN_NumEdgPerCrn = 4;
+  CrnNeigh~{1}~{1} = 2; CrnNeigh~{1}~{2} = 1; CrnNeigh~{1}~{3} = 3; CrnNeigh~{1}~{4} = 4;
+  CrnNeigh~{2}~{1} = 1; CrnNeigh~{2}~{2} = 5; CrnNeigh~{2}~{3} = 8; CrnNeigh~{2}~{4} = 9;
+  CrnNeigh~{3}~{1} = 2; CrnNeigh~{3}~{2} = 5; CrnNeigh~{3}~{3} = 6; CrnNeigh~{3}~{4} = 10;
+  CrnNeigh~{4}~{1} = 3; CrnNeigh~{4}~{2} = 6; CrnNeigh~{4}~{3} = 7; CrnNeigh~{4}~{4} = 11;
+  CrnNeigh~{5}~{1} = 4; CrnNeigh~{5}~{2} = 7; CrnNeigh~{5}~{3} = 8; CrnNeigh~{5}~{4} = 12;
+  CrnNeigh~{6}~{1} = 10; CrnNeigh~{6}~{2} = 9; CrnNeigh~{6}~{3} = 11; CrnNeigh~{6}~{4} = 12;
+EndIf
+
+If(FAC_Num == 12)
+  EDG_Num = 30;
+  EdgNeigh~{1}~{1} = 1; EdgNeigh~{1}~{2} = 2;
+  EdgNeigh~{2}~{1} = 1; EdgNeigh~{2}~{2} = 3;
+  EdgNeigh~{3}~{1} = 1; EdgNeigh~{3}~{2} = 4;
+  EdgNeigh~{4}~{1} = 1; EdgNeigh~{4}~{2} = 5;
+  EdgNeigh~{5}~{1} = 1; EdgNeigh~{5}~{2} = 6;
+  EdgNeigh~{6}~{1} = 2; EdgNeigh~{6}~{2} = 6;
+  EdgNeigh~{7}~{1} = 2; EdgNeigh~{7}~{2} = 3;
+  EdgNeigh~{8}~{1} = 3; EdgNeigh~{8}~{2} = 4;
+  EdgNeigh~{9}~{1} = 4; EdgNeigh~{9}~{2} = 5;
+  EdgNeigh~{10}~{1} = 5; EdgNeigh~{10}~{2} = 6;
+  EdgNeigh~{11}~{1} = 6; EdgNeigh~{11}~{2} = 7;
+  EdgNeigh~{12}~{1} = 2; EdgNeigh~{12}~{2} = 7;
+  EdgNeigh~{13}~{1} = 2; EdgNeigh~{13}~{2} = 8;
+  EdgNeigh~{14}~{1} = 3; EdgNeigh~{14}~{2} = 8;
+  EdgNeigh~{15}~{1} = 3; EdgNeigh~{15}~{2} = 9;
+  EdgNeigh~{16}~{1} = 4; EdgNeigh~{16}~{2} = 9;
+  EdgNeigh~{17}~{1} = 4; EdgNeigh~{17}~{2} = 10;
+  EdgNeigh~{18}~{1} = 5; EdgNeigh~{18}~{2} = 10;
+  EdgNeigh~{19}~{1} = 5; EdgNeigh~{19}~{2} = 11;
+  EdgNeigh~{20}~{1} = 6; EdgNeigh~{20}~{2} = 11;
+  EdgNeigh~{21}~{1} = 7; EdgNeigh~{21}~{2} = 8;
+  EdgNeigh~{22}~{1} = 8; EdgNeigh~{22}~{2} = 9;
+  EdgNeigh~{23}~{1} = 9; EdgNeigh~{23}~{2} = 10;
+  EdgNeigh~{24}~{1} = 10; EdgNeigh~{24}~{2} = 11;
+  EdgNeigh~{25}~{1} = 7; EdgNeigh~{25}~{2} = 11;
+  EdgNeigh~{26}~{1} = 8; EdgNeigh~{26}~{2} = 12;
+  EdgNeigh~{27}~{1} = 9; EdgNeigh~{27}~{2} = 12;
+  EdgNeigh~{28}~{1} = 10; EdgNeigh~{28}~{2} = 12;
+  EdgNeigh~{29}~{1} = 11; EdgNeigh~{29}~{2} = 12;
+  EdgNeigh~{30}~{1} = 7; EdgNeigh~{30}~{2} = 12;
+  CRN_Num = 20;
+  CRN_NumEdgPerCrn = 3;
+  CrnNeigh~{1}~{1}  =  1; CrnNeigh~{1}~{2}  =  5; CrnNeigh~{1}~{3}  =  6;
+  CrnNeigh~{2}~{1}  =  1; CrnNeigh~{2}~{2}  =  2; CrnNeigh~{2}~{3}  =  7;
+  CrnNeigh~{3}~{1}  =  2; CrnNeigh~{3}~{2}  =  3; CrnNeigh~{3}~{3}  =  8;
+  CrnNeigh~{4}~{1}  =  3; CrnNeigh~{4}~{2}  =  4; CrnNeigh~{4}~{3}  =  9;
+  CrnNeigh~{5}~{1}  =  4; CrnNeigh~{5}~{2}  =  5; CrnNeigh~{5}~{3}  = 10;
+  CrnNeigh~{6}~{1}  =  6; CrnNeigh~{6}~{2}  = 12; CrnNeigh~{6}~{3}  = 11;
+  CrnNeigh~{7}~{1}  =  7; CrnNeigh~{7}~{2}  = 13; CrnNeigh~{7}~{3}  = 14;
+  CrnNeigh~{8}~{1}  =  8; CrnNeigh~{8}~{2}  = 15; CrnNeigh~{8}~{3}  = 16;
+  CrnNeigh~{9}~{1}  =  9; CrnNeigh~{9}~{2}  = 17; CrnNeigh~{9}~{3}  = 18;
+  CrnNeigh~{10}~{1} = 10; CrnNeigh~{10}~{2} = 19; CrnNeigh~{10}~{3} = 20;
+  CrnNeigh~{11}~{1} = 12; CrnNeigh~{11}~{2} = 13; CrnNeigh~{11}~{3} = 21;
+  CrnNeigh~{12}~{1} = 14; CrnNeigh~{12}~{2} = 15; CrnNeigh~{12}~{3} = 22;
+  CrnNeigh~{13}~{1} = 16; CrnNeigh~{13}~{2} = 17; CrnNeigh~{13}~{3} = 23;
+  CrnNeigh~{14}~{1} = 18; CrnNeigh~{14}~{2} = 19; CrnNeigh~{14}~{3} = 24;
+  CrnNeigh~{15}~{1} = 11; CrnNeigh~{15}~{2} = 20; CrnNeigh~{15}~{3} = 25;
+  CrnNeigh~{16}~{1} = 21; CrnNeigh~{16}~{2} = 30; CrnNeigh~{16}~{3} = 26;
+  CrnNeigh~{17}~{1} = 22; CrnNeigh~{17}~{2} = 26; CrnNeigh~{17}~{3} = 27;
+  CrnNeigh~{18}~{1} = 23; CrnNeigh~{18}~{2} = 27; CrnNeigh~{18}~{3} = 28;
+  CrnNeigh~{19}~{1} = 24; CrnNeigh~{19}~{2} = 28; CrnNeigh~{19}~{3} = 29;
+  CrnNeigh~{20}~{1} = 25; CrnNeigh~{20}~{2} = 30; CrnNeigh~{20}~{3} = 29;
+EndIf
+
+If(FAC_Num == 20)
+  EDG_Num = 30;
+  EdgNeigh~{1}~{1} = 1; EdgNeigh~{1}~{2} = 5;
+  EdgNeigh~{2}~{1} = 1; EdgNeigh~{2}~{2} = 2;
+  EdgNeigh~{3}~{1} = 2; EdgNeigh~{3}~{2} = 3;
+  EdgNeigh~{4}~{1} = 3; EdgNeigh~{4}~{2} = 4;
+  EdgNeigh~{5}~{1} = 4; EdgNeigh~{5}~{2} = 5;
+  EdgNeigh~{6}~{1} = 1; EdgNeigh~{6}~{2} = 6;
+  EdgNeigh~{7}~{1} = 2; EdgNeigh~{7}~{2} = 7;
+  EdgNeigh~{8}~{1} = 3; EdgNeigh~{8}~{2} = 8;
+  EdgNeigh~{9}~{1} = 4; EdgNeigh~{9}~{2} = 9;
+  EdgNeigh~{10}~{1} = 5; EdgNeigh~{10}~{2} = 10;
+  EdgNeigh~{11}~{1} = 6; EdgNeigh~{11}~{2} = 15;
+  EdgNeigh~{12}~{1} = 6; EdgNeigh~{12}~{2} = 11;
+  EdgNeigh~{13}~{1} = 7; EdgNeigh~{13}~{2} = 11;
+  EdgNeigh~{14}~{1} = 7; EdgNeigh~{14}~{2} = 12;
+  EdgNeigh~{15}~{1} = 8; EdgNeigh~{15}~{2} = 12;
+  EdgNeigh~{16}~{1} = 8; EdgNeigh~{16}~{2} = 13;
+  EdgNeigh~{17}~{1} = 9; EdgNeigh~{17}~{2} = 13;
+  EdgNeigh~{18}~{1} = 9; EdgNeigh~{18}~{2} = 14;
+  EdgNeigh~{19}~{1} = 10; EdgNeigh~{19}~{2} = 14;
+  EdgNeigh~{20}~{1} = 10; EdgNeigh~{20}~{2} = 15;
+  EdgNeigh~{21}~{1} = 11; EdgNeigh~{21}~{2} = 16;
+  EdgNeigh~{22}~{1} = 12; EdgNeigh~{22}~{2} = 17;
+  EdgNeigh~{23}~{1} = 13; EdgNeigh~{23}~{2} = 18;
+  EdgNeigh~{24}~{1} = 14; EdgNeigh~{24}~{2} = 19;
+  EdgNeigh~{25}~{1} = 15; EdgNeigh~{25}~{2} = 20;
+  EdgNeigh~{26}~{1} = 16; EdgNeigh~{26}~{2} = 20;
+  EdgNeigh~{27}~{1} = 16; EdgNeigh~{27}~{2} = 17;
+  EdgNeigh~{28}~{1} = 17; EdgNeigh~{28}~{2} = 18;
+  EdgNeigh~{29}~{1} = 18; EdgNeigh~{29}~{2} = 19;
+  EdgNeigh~{30}~{1} = 19; EdgNeigh~{30}~{2} = 20;
+  CRN_Num = 12;
+  CRN_NumEdgPerCrn = 5;
+  CrnNeigh~{1}~{1}  =  1; CrnNeigh~{1}~{2}  =  2; CrnNeigh~{1}~{3}  =  3; CrnNeigh~{1}~{4}  =  4; CrnNeigh~{1}~{5}  =  5;
+  CrnNeigh~{2}~{1}  =  1; CrnNeigh~{2}~{2}  =  6; CrnNeigh~{2}~{3}  = 10; CrnNeigh~{2}~{4}  = 11; CrnNeigh~{2}~{5}  = 20;
+  CrnNeigh~{3}~{1}  =  2; CrnNeigh~{3}~{2}  =  6; CrnNeigh~{3}~{3}  =  7; CrnNeigh~{3}~{4}  = 12; CrnNeigh~{3}~{5}  = 13;
+  CrnNeigh~{4}~{1}  =  3; CrnNeigh~{4}~{2}  =  7; CrnNeigh~{4}~{3}  =  8; CrnNeigh~{4}~{4}  = 14; CrnNeigh~{4}~{5}  = 15;
+  CrnNeigh~{5}~{1}  =  4; CrnNeigh~{5}~{2}  =  8; CrnNeigh~{5}~{3}  =  9; CrnNeigh~{5}~{4}  = 16; CrnNeigh~{5}~{5}  = 17;
+  CrnNeigh~{6}~{1}  =  5; CrnNeigh~{6}~{2}  =  9; CrnNeigh~{6}~{3}  = 10; CrnNeigh~{6}~{4}  = 18; CrnNeigh~{6}~{5}  = 19;
+  CrnNeigh~{7}~{1}  = 11; CrnNeigh~{7}~{2}  = 12; CrnNeigh~{7}~{3}  = 21; CrnNeigh~{7}~{4}  = 25; CrnNeigh~{7}~{5}  = 26;
+  CrnNeigh~{8}~{1}  = 13; CrnNeigh~{8}~{2}  = 14; CrnNeigh~{8}~{3}  = 21; CrnNeigh~{8}~{4}  = 22; CrnNeigh~{8}~{5}  = 27;
+  CrnNeigh~{9}~{1}  = 15; CrnNeigh~{9}~{2}  = 16; CrnNeigh~{9}~{3}  = 22; CrnNeigh~{9}~{4}  = 23; CrnNeigh~{9}~{5}  = 28;
+  CrnNeigh~{10}~{1} = 17; CrnNeigh~{10}~{2} = 18; CrnNeigh~{10}~{3} = 23; CrnNeigh~{10}~{4} = 24; CrnNeigh~{10}~{5} = 29;
+  CrnNeigh~{11}~{1} = 19; CrnNeigh~{11}~{2} = 20; CrnNeigh~{11}~{3} = 24; CrnNeigh~{11}~{4} = 25; CrnNeigh~{11}~{5} = 30;
+  CrnNeigh~{12}~{1} = 26; CrnNeigh~{12}~{2} = 27; CrnNeigh~{12}~{3} = 28; CrnNeigh~{12}~{4} = 29; CrnNeigh~{12}~{5} = 30;
+EndIf
+
+If(FAC_Num == 1)
+  EDG_Num = 0;
+  EDG_NumFacPerEdg = 0;
+  CRN_Num = 0;
+  CRN_NumEdgPerCrn = 0;
+Else
+  EDG_NumFacPerEdg = 2;
+EndIf
+
+For i In {1:CRN_Num}
+  CRN~{i} = i;
+EndFor
+For i In {1:EDG_Num}
+  EDG~{i} = 100+i;
+EndFor
+For i In {1:FAC_Num}
+  SUR~{i} = 200+i;
+EndFor
+SUR_Scatt = 200;
+VOL = 301;

HelmholtzHABCwithCorners/padePolyhedron.geo

+Include "padePolyhedron.dat";
+
+SetFactory("OpenCASCADE");
+
+/// TRUNCATED POLYHEDRAL DOMAIN
+
+X_SCA = 0.;
+Y_SCA = 0.;
+Z_SCA = 0.;
+
+If(FAC_Num == 4)
+  dimL = dimRext/Sqrt[3];
+
+  P_1 = newp; Point(P_1) = { dimL-X_SCA, dimL-Y_SCA, dimL-Z_SCA};
+  P_2 = newp; Point(P_2) = { dimL-X_SCA,-dimL-Y_SCA,-dimL-Z_SCA};
+  P_3 = newp; Point(P_3) = {-dimL-X_SCA, dimL-Y_SCA,-dimL-Z_SCA};
+  P_4 = newp; Point(P_4) = {-dimL-X_SCA,-dimL-Y_SCA, dimL-Z_SCA};
+
+  L_1 = newl; Line(L_1) = {P_1, P_2};
+  L_2 = newl; Line(L_2) = {P_1, P_3};
+  L_3 = newl; Line(L_3) = {P_1, P_4};
+  L_4 = newl; Line(L_4) = {P_2, P_3};
+  L_5 = newl; Line(L_5) = {P_2, P_4};
+  L_6 = newl; Line(L_6) = {P_3, P_4};
+
+  LL_1 = newll; Line Loop(LL_1) = {L_1, L_2, L_4};
+  LL_2 = newll; Line Loop(LL_2) = {L_3, L_1, L_5};
+  LL_3 = newll; Line Loop(LL_3) = {L_2, L_3, L_6};
+  LL_4 = newll; Line Loop(LL_4) = {L_5, L_4, L_6};
+EndIf
+
+If(FAC_Num == 6)
+  dimL = dimRext/Sqrt[3];
+
+  P_1 = newp; Point(P_1) = {-dimL-X_SCA,-dimL-Y_SCA,-dimL-Z_SCA};
+  P_2 = newp; Point(P_2) = { dimL-X_SCA,-dimL-Y_SCA,-dimL-Z_SCA};
+  P_3 = newp; Point(P_3) = { dimL-X_SCA,-dimL-Y_SCA, dimL-Z_SCA};
+  P_4 = newp; Point(P_4) = {-dimL-X_SCA,-dimL-Y_SCA, dimL-Z_SCA};
+  P_5 = newp; Point(P_5) = {-dimL-X_SCA, dimL-Y_SCA, dimL-Z_SCA};
+  P_6 = newp; Point(P_6) = { dimL-X_SCA, dimL-Y_SCA, dimL-Z_SCA};
+  P_7 = newp; Point(P_7) = {-dimL-X_SCA, dimL-Y_SCA,-dimL-Z_SCA};
+  P_8 = newp; Point(P_8) = { dimL-X_SCA, dimL-Y_SCA,-dimL-Z_SCA};
+
+  L_1 = newl; Line(L_1) = {P_1, P_2};
+  L_2 = newl; Line(L_2) = {P_4, P_3};
+  L_3 = newl; Line(L_3) = {P_5, P_6};
+  L_4 = newl; Line(L_4) = {P_7, P_8};
+  L_5 = newl; Line(L_5) = {P_1, P_7};
+  L_6 = newl; Line(L_6) = {P_2, P_8};
+  L_7 = newl; Line(L_7) = {P_3, P_6};
+  L_8 = newl; Line(L_8) = {P_4, P_5};
+  L_9 = newl; Line(L_9) = {P_1, P_4};
+  L_10 = newl; Line(L_10) = {P_7, P_5};
+  L_11 = newl; Line(L_11) = {P_8, P_6};
+  L_12 = newl; Line(L_12) = {P_2, P_3};
+
+  LL_1 = newll; Line Loop(LL_1) = {L_9, L_8, L_10, L_5};
+  LL_2 = newll; Line Loop(LL_2) = {L_6, L_11, L_7, L_12};
+  LL_3 = newll; Line Loop(LL_3) = {L_1, L_12, L_2, L_9};
+  LL_4 = newll; Line Loop(LL_4) = {L_10, L_3, L_11, L_4};
+  LL_5 = newll; Line Loop(LL_5) = {L_5, L_4, L_6, L_1};
+  LL_6 = newll; Line Loop(LL_6) = {L_2, L_7, L_3, L_8};
+EndIf
+
+If(FAC_Num == 8)
+  dimL = dimRext;
+
+  P_1 = newp; Point(P_1) = {-dimL-X_SCA,     -Y_SCA,     -Z_SCA};
+  P_2 = newp; Point(P_2) = {     -X_SCA,-dimL-Y_SCA,     -Z_SCA};
+  P_3 = newp; Point(P_3) = {     -X_SCA,     -Y_SCA,-dimL-Z_SCA};
+  P_4 = newp; Point(P_4) = {     -X_SCA, dimL-Y_SCA,     -Z_SCA};
+  P_5 = newp; Point(P_5) = {     -X_SCA,     -Y_SCA, dimL-Z_SCA};
+  P_6 = newp; Point(P_6) = { dimL-X_SCA,     -Y_SCA,     -Z_SCA};
+
+  L_1 = newl; Line(L_1) = {P_2, P_1};
+  L_2 = newl; Line(L_2) = {P_3, P_1};
+  L_3 = newl; Line(L_3) = {P_4, P_1};
+  L_4 = newl; Line(L_4) = {P_5, P_1};
+  L_5 = newl; Line(L_5) = {P_2, P_3};
+  L_6 = newl; Line(L_6) = {P_3, P_4};
+  L_7 = newl; Line(L_7) = {P_4, P_5};
+  L_8 = newl; Line(L_8) = {P_5, P_2};
+  L_9 = newl; Line(L_9) = {P_6, P_2};
+  L_10 = newl; Line(L_10) = {P_6, P_3};
+  L_11 = newl; Line(L_11) = {P_6, P_4};
+  L_12 = newl; Line(L_12) = {P_6, P_5};
+
+  LL_1 = newll; Line Loop(LL_1) = {L_1, L_2, L_5};
+  LL_2 = newll; Line Loop(LL_2) = {L_2, L_3, L_6};
+  LL_3 = newll; Line Loop(LL_3) = {L_3, L_4, L_7};
+  LL_4 = newll; Line Loop(LL_4) = {L_4, L_1, L_8};
+  LL_5 = newll; Line Loop(LL_5) = {L_9, L_10, L_5};
+  LL_6 = newll; Line Loop(LL_6) = {L_10, L_11, L_6};
+  LL_7 = newll; Line Loop(LL_7) = {L_11, L_12, L_7};
+  LL_8 = newll; Line Loop(LL_8) = {L_12, L_9, L_8};
+EndIf
+
+If(FAC_Num == 12)
+  dimL = dimRext/Tan[Pi/3];
+  GoRat = 1.61803398875*dimL;
+  invGR = 0.61803398875*dimL;
+
+  P_1  = newp; Point(P_1)  = { -dimL-X_SCA, -dimL-Y_SCA, -dimL-Z_SCA};
+  P_2  = newp; Point(P_2)  = {-invGR-X_SCA,      -Y_SCA,-GoRat-Z_SCA};
+  P_3  = newp; Point(P_3)  = { invGR-X_SCA,      -Y_SCA,-GoRat-Z_SCA};
+  P_4  = newp; Point(P_4)  = {  dimL-X_SCA, -dimL-Y_SCA, -dimL-Z_SCA};
+  P_5  = newp; Point(P_5)  = {      -X_SCA,-GoRat-Y_SCA,-invGR-Z_SCA};
+  P_6  = newp; Point(P_6)  = {-GoRat-X_SCA,-invGR-Y_SCA,      -Z_SCA};
+  P_7  = newp; Point(P_7)  = { -dimL-X_SCA,  dimL-Y_SCA, -dimL-Z_SCA};
+  P_8  = newp; Point(P_8)  = {  dimL-X_SCA,  dimL-Y_SCA, -dimL-Z_SCA};
+  P_9  = newp; Point(P_9)  = { GoRat-X_SCA,-invGR-Y_SCA,      -Z_SCA};
+  P_10 = newp; Point(P_10) = {      -X_SCA,-GoRat-Y_SCA, invGR-Z_SCA};
+  P_11 = newp; Point(P_11) = {-GoRat-X_SCA, invGR-Y_SCA,      -Z_SCA};
+  P_12 = newp; Point(P_12) = {      -X_SCA, GoRat-Y_SCA,-invGR-Z_SCA};
+  P_13 = newp; Point(P_13) = { GoRat-X_SCA, invGR-Y_SCA,      -Z_SCA};
+  P_14 = newp; Point(P_14) = {  dimL-X_SCA, -dimL-Y_SCA,  dimL-Z_SCA};
+  P_15 = newp; Point(P_15) = { -dimL-X_SCA, -dimL-Y_SCA,  dimL-Z_SCA};
+  P_16 = newp; Point(P_16) = { -dimL-X_SCA,  dimL-Y_SCA,  dimL-Z_SCA};
+  P_17 = newp; Point(P_17) = {      -X_SCA, GoRat-Y_SCA, invGR-Z_SCA};
+  P_18 = newp; Point(P_18) = {  dimL-X_SCA,  dimL-Y_SCA,  dimL-Z_SCA};
+  P_19 = newp; Point(P_19) = { invGR-X_SCA,      -Y_SCA, GoRat-Z_SCA};
+  P_20 = newp; Point(P_20) = {-invGR-X_SCA,      -Y_SCA, GoRat-Z_SCA};
+
+  L_1 = newl;  Line(L_1)  = {P_1, P_2};
+  L_2 = newl;  Line(L_2)  = {P_2, P_3};
+  L_3 = newl;  Line(L_3)  = {P_3, P_4};
+  L_4 = newl;  Line(L_4)  = {P_4, P_5};
+  L_5 = newl;  Line(L_5)  = {P_5, P_1};
+  L_6 = newl;  Line(L_6)  = {P_1, P_6};
+  L_7 = newl;  Line(L_7)  = {P_2, P_7};
+  L_8 = newl;  Line(L_8)  = {P_3, P_8};
+  L_9 = newl;  Line(L_9)  = {P_4, P_9};
+  L_10 = newl; Line(L_10) = {P_5, P_10};
+  L_11 = newl; Line(L_11) = {P_6, P_15};
+  L_12 = newl; Line(L_12) = {P_6, P_11};
+  L_13 = newl; Line(L_13) = {P_7, P_11};
+  L_14 = newl; Line(L_14) = {P_7, P_12};
+  L_15 = newl; Line(L_15) = {P_8, P_12};
+  L_16 = newl; Line(L_16) = {P_8, P_13};
+  L_17 = newl; Line(L_17) = {P_9, P_13};
+  L_18 = newl; Line(L_18) = {P_9, P_14};
+  L_19 = newl; Line(L_19) = {P_10, P_14};
+  L_20 = newl; Line(L_20) = {P_10, P_15};
+  L_21 = newl; Line(L_21) = {P_11, P_16};
+  L_22 = newl; Line(L_22) = {P_12, P_17};
+  L_23 = newl; Line(L_23) = {P_13, P_18};
+  L_24 = newl; Line(L_24) = {P_14, P_19};
+  L_25 = newl; Line(L_25) = {P_15, P_20};
+  L_26 = newl; Line(L_26) = {P_17, P_16};
+  L_27 = newl; Line(L_27) = {P_18, P_17};
+  L_28 = newl; Line(L_28) = {P_19, P_18};
+  L_29 = newl; Line(L_29) = {P_20, P_19};
+  L_30 = newl; Line(L_30) = {P_16, P_20};
+
+  LL_1  = newll; Line Loop(LL_1 ) = {L_1, L_2, L_3, L_4, L_5};
+  LL_2  = newll; Line Loop(LL_2 ) = {L_6, L_1, L_7, L_12, L_13};
+  LL_3  = newll; Line Loop(LL_3 ) = {L_7, L_2, L_8, L_14, L_15};
+  LL_4  = newll; Line Loop(LL_4 ) = {L_8, L_3, L_9, L_16, L_17};
+  LL_5  = newll; Line Loop(LL_5 ) = {L_9, L_4, L_10, L_18, L_19};
+  LL_6  = newll; Line Loop(LL_6 ) = {L_10, L_5, L_6, L_20, L_11};
+  LL_7  = newll; Line Loop(LL_7 ) = {L_11, L_12, L_25, L_21, L_30};
+  LL_8  = newll; Line Loop(LL_8 ) = {L_13, L_14, L_21, L_22, L_26};
+  LL_9  = newll; Line Loop(LL_9 ) = {L_15, L_16, L_22, L_23, L_27};
+  LL_10 = newll; Line Loop(LL_10) = {L_17, L_18, L_23, L_24, L_28};
+  LL_11 = newll; Line Loop(LL_11) = {L_19, L_20, L_24, L_25, L_29};
+  LL_12 = newll; Line Loop(LL_12) = {L_26, L_27, L_28, L_29, L_30};
+EndIf
+
+If(FAC_Num == 20)
+  dimL = dimRext/(Tan[Pi/5]*2.61803398875);
+  GoRat = 1.61803398875*dimL;
+
+  P_1  = newp; Point(P_1)  = {-GoRat-X_SCA,      -Y_SCA, -dimL-Z_SCA};
+  P_2  = newp; Point(P_2)  = {      -X_SCA,  dimL-Y_SCA,-GoRat-Z_SCA};
+  P_3  = newp; Point(P_3)  = {      -X_SCA, -dimL-Y_SCA,-GoRat-Z_SCA};
+  P_4  = newp; Point(P_4)  = { -dimL-X_SCA,-GoRat-Y_SCA,      -Z_SCA};
+  P_5  = newp; Point(P_5)  = {-GoRat-X_SCA,      -Y_SCA,  dimL-Z_SCA};
+  P_6  = newp; Point(P_6)  = { -dimL-X_SCA, GoRat-Y_SCA,      -Z_SCA};
+  P_7  = newp; Point(P_7)  = { GoRat-X_SCA,      -Y_SCA, -dimL-Z_SCA};
+  P_8  = newp; Point(P_8)  = {  dimL-X_SCA,-GoRat-Y_SCA,      -Z_SCA};
+  P_9  = newp; Point(P_9)  = {      -X_SCA, -dimL-Y_SCA, GoRat-Z_SCA};
+  P_10 = newp; Point(P_10) = {      -X_SCA,  dimL-Y_SCA, GoRat-Z_SCA};
+  P_11 = newp; Point(P_11) = {  dimL-X_SCA, GoRat-Y_SCA,      -Z_SCA};
+  P_12 = newp; Point(P_12) = { GoRat-X_SCA,      -Y_SCA,  dimL-Z_SCA};
+
+  L_1  = newl; Line(L_1)  = {P_1, P_2};
+  L_2  = newl; Line(L_2)  = {P_1, P_3};
+  L_3  = newl; Line(L_3)  = {P_1, P_4};
+  L_4  = newl; Line(L_4)  = {P_1, P_5};
+  L_5  = newl; Line(L_5)  = {P_1, P_6};
+  L_6  = newl; Line(L_6)  = {P_2, P_3};
+  L_7  = newl; Line(L_7)  = {P_3, P_4};
+  L_8  = newl; Line(L_8)  = {P_4, P_5};
+  L_9  = newl; Line(L_9)  = {P_5, P_6};
+  L_10 = newl; Line(L_10) = {P_6, P_2};
+  L_11 = newl; Line(L_11) = {P_2, P_7};
+  L_12 = newl; Line(L_12) = {P_7, P_3};
+  L_13 = newl; Line(L_13) = {P_3, P_8};
+  L_14 = newl; Line(L_14) = {P_8, P_4};
+  L_15 = newl; Line(L_15) = {P_4, P_9};
+  L_16 = newl; Line(L_16) = {P_9, P_5};
+  L_17 = newl; Line(L_17) = {P_5, P_10};
+  L_18 = newl; Line(L_18) = {P_10, P_6};
+  L_19 = newl; Line(L_19) = {P_6, P_11};
+  L_20 = newl; Line(L_20) = {P_11, P_2};
+  L_21 = newl; Line(L_21) = {P_7, P_8};
+  L_22 = newl; Line(L_22) = {P_8, P_9};
+  L_23 = newl; Line(L_23) = {P_9, P_10};
+  L_24 = newl; Line(L_24) = {P_10, P_11};
+  L_25 = newl; Line(L_25) = {P_11, P_7};
+  L_26 = newl; Line(L_26) = {P_7, P_12};
+  L_27 = newl; Line(L_27) = {P_8, P_12};
+  L_28 = newl; Line(L_28) = {P_9, P_12};
+  L_29 = newl; Line(L_29) = {P_10, P_12};
+  L_30 = newl; Line(L_30) = {P_11, P_12};
+
+  LL_1  = newll; Line Loop(LL_1 ) = {L_1, L_2, L_6};
+  LL_2  = newll; Line Loop(LL_2 ) = {L_2, L_3, L_7};
+  LL_3  = newll; Line Loop(LL_3 ) = {L_3, L_4, L_8};
+  LL_4  = newll; Line Loop(LL_4 ) = {L_4, L_5, L_9};
+  LL_5  = newll; Line Loop(LL_5 ) = {L_5, L_1, L_10};
+  LL_6  = newll; Line Loop(LL_6 ) = {L_11, L_6, L_12};
+  LL_7  = newll; Line Loop(LL_7 ) = {L_13, L_7, L_14};
+  LL_8  = newll; Line Loop(LL_8 ) = {L_15, L_8, L_16};
+  LL_9  = newll; Line Loop(LL_9 ) = {L_17, L_9, L_18};
+  LL_10 = newll; Line Loop(LL_10) = {L_19, L_10, L_20};
+  LL_11 = newll; Line Loop(LL_11) = {L_21, L_12, L_13};
+  LL_12 = newll; Line Loop(LL_12) = {L_22, L_14, L_15};
+  LL_13 = newll; Line Loop(LL_13) = {L_23, L_16, L_17};
+  LL_14 = newll; Line Loop(LL_14) = {L_24, L_18, L_19};
+  LL_15 = newll; Line Loop(LL_15) = {L_25, L_20, L_11};
+  LL_16 = newll; Line Loop(LL_16) = {L_26, L_21, L_27};
+  LL_17 = newll; Line Loop(LL_17) = {L_27, L_22, L_28};
+  LL_18 = newll; Line Loop(LL_18) = {L_28, L_23, L_29};
+  LL_19 = newll; Line Loop(LL_19) = {L_29, L_24, L_30};
+  LL_20 = newll; Line Loop(LL_20) = {L_30, L_25, L_26};
+EndIf
+
+/// EXTERIOR BOUNDARY
+
+VOL_Ext = newv;
+If(FAC_Num == 1)
+  Sphere(VOL_Ext) = {0.,0.,0.,dimRmean};
+  Physical Surface(SUR~{1}) = { CombinedBoundary{ Volume{VOL_Ext}; }};
+Else
+  For i In {1:FAC_Num}
+    S~{i} = news; Plane Surface(S~{i}) = {LL~{i}};
+    listS[] += S~{i};
+  EndFor
+  SL = newsl; Surface Loop(SL) = {listS[]};
+  Volume(VOL_Ext) = {SL};
+
+  For i In {1:CRN_Num}
+    Physical Point(CRN~{i}) = {P~{i}};
+  EndFor
+  For i In {1:EDG_Num}
+    Physical Line(EDG~{i}) = {L~{i}};
+  EndFor
+  For i In {1:FAC_Num}
+    Physical Surface(SUR~{i}) = {S~{i}};
+  EndFor
+EndIf
+
+/// SCATTERER BOUNDARY
+
+VOL_Scatt = newv;
+Sphere(VOL_Scatt) = {0.,0.,0.,R_SCA};
+Physical Surface(SUR_Scatt) = { CombinedBoundary{ Volume{VOL_Scatt}; }};
+
+/// MAIN DOMAIN
+
+VOL_Dom = newv;
+BooleanDifference(VOL_Dom) = { Volume{VOL_Ext}; }{ Volume{VOL_Scatt}; };
+Physical Volume(VOL) = {VOL_Dom};
+Delete{ Volume{VOL_Ext,VOL_Scatt}; }
+
+/// GENERATE MESHES
+
+SetOrder 1;
+Mesh.ElementOrder = 1;
+Mesh 2;
+Save "mainCurv.msh";
+SetOrder ORDER;
+Mesh.ElementOrder = ORDER;

HelmholtzHABCwithCorners/padeSquare.dat

+BND_Neumann        = 0;
+BND_Sommerfeld     = 1;
+BND_Second         = 2;
+BND_Pade           = 3;
+BND_CRBC           = 4;
+BND_PML            = 5;
+
+DefineConstant[
+  BND_TYPE = {BND_Pade,
+    Name "Input/5Model/02Boundary condition (edges)", Highlight "Blue",
+    Choices {BND_Neumann    = "Homogeneous Neumann",
+             BND_Sommerfeld = "Sommerfeld ABC",
+             BND_Second     = "Second-order ABC",
+             BND_Pade       = "Pade ABC",
+             BND_CRBC       = "CRBC",
+             BND_PML        = "PML"}}
+];
+
+DefineConstant[
+  dimL  = { 2.2, Min 1, Step 0.1,  Max 10, Name "Input/1Geometry/1Domain length"},
+  X_SCA = { 1.1, Min 1, Step 0.01, Max 10, Name "Input/1Geometry/2Scatterer position (x0)"},
+  Y_SCA = { 1.1, Min 1, Step 0.01, Max 10, Name "Input/1Geometry/3Scatterer position (y0)"},
+  Npml  = { 3,   Min 1, Step 1,    Max 5, Name "Input/5Model/03PML: Thickness (N*Lc)", Visible (BND_TYPE == BND_PML)},
+  //sigmaPml  = { 2e4, Choices {2e4}, Name "Input/5Model/03PML: Sigma Mult", Visible (BND_TYPE == BND_PML)},
+  rotPml  = { 0, Min -1.57079632679, Max 1.57079632679, Step 0.0001, Name "Input/5Model/03PML: Rotation", Visible (BND_TYPE == BND_PML)}
+]; // 1.57079632679
+
+Lpml = LC*Npml;
+
+X~{1} = -X_SCA;
+Y~{1} = -Y_SCA;
+X~{2} = -X_SCA+dimL;
+Y~{2} = -Y_SCA;
+X~{3} = -X_SCA+dimL;
+Y~{3} = -Y_SCA+dimL;
+X~{4} = -X_SCA;
+Y~{4} = -Y_SCA+dimL;
+
+CRN_1 = 101;
+CRN_2 = 102;
+CRN_3 = 103;
+CRN_4 = 104;
+
+BND_1 = 201;
+BND_2 = 202;
+BND_3 = 203;
+BND_4 = 204;
+
+BND_Scatt = 205;
+BND_PmlExt = 206;
+
+DOM = 301;
+DOM_PML = 302;