diff --git a/HelmholtzHABCwithCorners/README.txt b/HelmholtzHABCwithCorners/README.txt
new file mode 100644
index 0000000000000000000000000000000000000000..71ad918d0d76348621fc39a002641dc3cbbeccb1
--- /dev/null
+++ b/HelmholtzHABCwithCorners/README.txt
@@ -0,0 +1,11 @@
+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
diff --git a/HelmholtzHABCwithCorners/main.dat b/HelmholtzHABCwithCorners/main.dat
new file mode 100644
index 0000000000000000000000000000000000000000..27607b914940198ba734d7cdaaf59fbcddb6ac95
--- /dev/null
+++ b/HelmholtzHABCwithCorners/main.dat
@@ -0,0 +1,77 @@
+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
diff --git a/HelmholtzHABCwithCorners/main.geo b/HelmholtzHABCwithCorners/main.geo
new file mode 100644
index 0000000000000000000000000000000000000000..ace41bb18e373da98c43ac619a559bc33831f825
--- /dev/null
+++ b/HelmholtzHABCwithCorners/main.geo
@@ -0,0 +1,15 @@
+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];
diff --git a/HelmholtzHABCwithCorners/main.pro b/HelmholtzHABCwithCorners/main.pro
new file mode 100644
index 0000000000000000000000000000000000000000..d9329fbd0bac321b9f3252ba8d47134d4a4aab34
--- /dev/null
+++ b/HelmholtzHABCwithCorners/main.pro
@@ -0,0 +1,122 @@
+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 }
+];
diff --git a/HelmholtzHABCwithCorners/padeCube.dat b/HelmholtzHABCwithCorners/padeCube.dat
new file mode 100644
index 0000000000000000000000000000000000000000..dede378c221ad8e06879f9ca60d668fa2bdf6bb1
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padeCube.dat
@@ -0,0 +1,36 @@
+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;
diff --git a/HelmholtzHABCwithCorners/padeCube.geo b/HelmholtzHABCwithCorners/padeCube.geo
new file mode 100644
index 0000000000000000000000000000000000000000..7e470a40b086052d1503d45d5b41556254b8ee67
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padeCube.geo
@@ -0,0 +1,121 @@
+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};
diff --git a/HelmholtzHABCwithCorners/padeCube.pro b/HelmholtzHABCwithCorners/padeCube.pro
new file mode 100644
index 0000000000000000000000000000000000000000..76d8165f7b257f73bfd1ebbd369cff57dd6f9e7d
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padeCube.pro
@@ -0,0 +1,508 @@
+Include "padeCube.dat";
+
+//==================================================================================================
+// OPTIONS and PARAMETERS
+//==================================================================================================
+
+BND_Neumann = 0;
+BND_Sommerfeld = 1;
+BND_Second = 2;
+BND_Pade = 3;
+CRN_Regularization = 0;
+CRN_Compatibility = 1;
+KEPS_Nothing = 0;
+KEPS_Analytic = 1;
+KEPS_Numeric = 2;
+
+DefineConstant[
+ BND_TYPE = {BND_Pade,
+ Name "Input/5Model/02Boundary condition (faces)", Highlight "Blue",
+ Choices {BND_Neumann = "Homogeneous Neumann",
+ BND_Sommerfeld = "Sommerfeld ABC",
+ BND_Second = "Second-order ABC",
+ BND_Pade = "Pade ABC"}},
+ CRN_TYPE = {CRN_Compatibility,
+ Name "Input/5Model/03Boundary condition (edges, corners)", Highlight "Red",
+ Visible ((BND_TYPE == BND_Second) || (BND_TYPE == BND_Pade)),
+ Choices {CRN_Regularization = "Regularization",
+ CRN_Compatibility = "Compatibility"}},
+ nPade = {4, Min 0, Step 1, Max 6,
+ Name "Input/5Model/05Pade: Number of fields",
+ Visible (BND_TYPE == BND_Pade)},
+ thetaPadeInput = {3, Min 0, Step 1, Max 4,
+ Name "Input/5Model/06Pade: Rotation of branch cut",
+ Visible (BND_TYPE == BND_Pade)},
+ KEPS_TYPE = {KEPS_Nothing,
+ Name "Input/5Model/07Curvature for regularization", Highlight "Red",
+ Visible ((BND_TYPE == BND_Pade) && (CRN_TYPE == CRN_Regularization)),
+ Choices {KEPS_Nothing = "Nothing",
+ KEPS_Analytic = "Analytic formula",
+ KEPS_Numeric = "Numerical curvature"}}
+];
+
+If(BND_TYPE == BND_Pade)
+ If(thetaPadeInput == 0)
+ thetaPade = 0;
+ EndIf
+ If(thetaPadeInput == 1)
+ thetaPade = Pi/8;
+ EndIf
+ If(thetaPadeInput == 2)
+ thetaPade = Pi/4;
+ EndIf
+ If(thetaPadeInput == 3)
+ thetaPade = Pi/3;
+ EndIf
+ If(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
+
+Group {
+ Dom = Region[{VOL}];
+ BndSca = Region[{SUR_Scatt}];
+
+ For i In {1:2}
+ FacesAll += Region[{SUR~{i}~{0}~{0}}];
+ FacesAll += Region[{SUR~{0}~{i}~{0}}];
+ FacesAll += Region[{SUR~{0}~{0}~{i}}];
+ FacesX += Region[{SUR~{i}~{0}~{0}}];
+ FacesY += Region[{SUR~{0}~{i}~{0}}];
+ FacesZ += Region[{SUR~{0}~{0}~{i}}];
+ EndFor
+ For i In {1:2}
+ For j In {1:2}
+ EdgesAll += Region[{LIN~{0}~{i}~{j}}];
+ EdgesAll += Region[{LIN~{j}~{0}~{i}}];
+ EdgesAll += Region[{LIN~{i}~{j}~{0}}];
+ EdgesYZ += Region[{LIN~{0}~{i}~{j}}];
+ EdgesZX += Region[{LIN~{j}~{0}~{i}}];
+ EdgesXY += Region[{LIN~{i}~{j}~{0}}];
+ BndFacesX += Region[{LIN~{j}~{0}~{i}}];
+ BndFacesX += Region[{LIN~{i}~{j}~{0}}];
+ BndFacesY += Region[{LIN~{0}~{i}~{j}}];
+ BndFacesY += Region[{LIN~{i}~{j}~{0}}];
+ BndFacesZ += Region[{LIN~{0}~{i}~{j}}];
+ BndFacesZ += Region[{LIN~{j}~{0}~{i}}];
+ EndFor
+ EndFor
+ For i In {1:2}
+ For j In {1:2}
+ For k In {1:2}
+ CornersAll += Region[{PNT~{i}~{j}~{k}}];
+ BndEdgesYZ += Region[{PNT~{i}~{j}~{k}}];
+ BndEdgesZX += Region[{PNT~{i}~{j}~{k}}];
+ BndEdgesXY += Region[{PNT~{i}~{j}~{k}}];
+ EndFor
+ EndFor
+ EndFor
+
+ BndExt = Region[{FacesAll,EdgesAll,CornersAll}];
+ DomAll = Region[{Dom,BndSca,BndExt}];
+}
+
+Function {
+ NormalNum[] = VectorField[XYZ[]]{1001};
+ CurvNum[] = ScalarField[XYZ[]]{1002};
+
+If((BND_TYPE == BND_Pade) && (CRN_TYPE == CRN_Regularization) && (KEPS_TYPE == KEPS_Numeric))
+ kEps[BndExt] = WAVENUMBER + I[] * 0.39 * WAVENUMBER^(1/3) * (CurvNum[])^(2/3);
+Else
+ kEps[BndExt] = WAVENUMBER;
+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.]];
+EndIf
+}
+
+//==================================================================================================
+// FONCTION SPACES with CONSTRAINTS
+//==================================================================================================
+
+Constraint {
+ { Name DirichletBC;
+ Case {{ Region BndSca; Value f_ref[]; }}
+ }
+}
+
+FunctionSpace {
+ { Name H_nx; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{FacesAll,BndSca}]; Entity NodesOf[All]; }}}
+ { Name H_ny; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{FacesAll,BndSca}]; Entity NodesOf[All]; }}}
+ { Name H_nz; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{FacesAll,BndSca}]; Entity NodesOf[All]; }}}
+ { Name H_cur; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{FacesAll,BndSca}]; Entity NodesOf[All]; }}}
+ { Name H_num; Type Form0;
+ BasisFunction {{ Name sn; NameOfCoef pn; Function BF_Node; Support Region[{DomAll}]; Entity NodesOf[All]; }}
+ Constraint {{ NameOfCoef pn; EntityType NodesOf; NameOfConstraint DirichletBC; }}
+ }
+If(BND_TYPE == BND_Pade)
+If(CRN_TYPE == CRN_Regularization)
+ For m In {1:nPade}
+ { Name H~{m}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{FacesAll}]; Entity NodesOf[All]; }}}
+ EndFor
+EndIf
+If(CRN_TYPE == CRN_Compatibility)
+ For m In {1:nPade}
+ { Name H~{m}~{0}~{0}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{FacesX,BndFacesX}]; Entity NodesOf[All]; }}}
+ { Name H~{0}~{m}~{0}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{FacesY,BndFacesY}]; Entity NodesOf[All]; }}}
+ { Name H~{0}~{0}~{m}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{FacesZ,BndFacesZ}]; Entity NodesOf[All]; }}}
+ EndFor
+ For m In {1:nPade}
+ For n In {1:nPade}
+ { Name H~{0}~{m}~{n}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{EdgesYZ,BndEdgesYZ}]; Entity NodesOf[All]; }}}
+ { Name H~{n}~{0}~{m}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{EdgesZX,BndEdgesZX}]; Entity NodesOf[All]; }}}
+ { Name H~{m}~{n}~{0}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{EdgesXY,BndEdgesXY}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+ For m In {1:nPade}
+ For n In {1:nPade}
+ For o In {1:nPade}
+ { Name H~{m}~{n}~{o}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{CornersAll}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+ EndFor
+EndIf
+EndIf
+}
+
+//==================================================================================================
+// FORMULATIONS
+//==================================================================================================
+
+Formulation {
+ { Name NumNormal; Type FemEquation;
+ Quantity {
+ { Name u_nx; Type Local; NameOfSpace H_nx; }
+ { Name u_ny; Type Local; NameOfSpace H_ny; }
+ { Name u_nz; Type Local; NameOfSpace H_nz; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_nx} , {u_nx} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ny} , {u_ny} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_nz} , {u_nz} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompX[Normal[]] , {u_nx} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompY[Normal[]] , {u_ny} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompZ[Normal[]] , {u_nz} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ }
+ }
+
+ { Name NumCur; Type FemEquation;
+ Quantity {
+ { Name u_nx; Type Local; NameOfSpace H_nx; }
+ { Name u_ny; Type Local; NameOfSpace H_ny; }
+ { Name u_nz; Type Local; NameOfSpace H_nz; }
+ { Name u_cur; Type Local; NameOfSpace H_cur; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_nx} , {u_nx} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ny} , {u_ny} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_nz} , {u_nz} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompX[NormalNum[]] , {u_nx} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompY[NormalNum[]] , {u_ny} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompZ[NormalNum[]] , {u_nz} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+
+ Galerkin{ [ Dof{u_cur} , {u_cur} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[0.5,0,0] * Dof{d u_nx} , {u_cur} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[0,0.5,0] * Dof{d u_ny} , {u_cur} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[0,0,0.5] * Dof{d u_nz} , {u_cur} ]; In Region[{BndExt,BndSca}]; Jacobian JSur; Integration I1; }
+ }
+ }
+
+ { Name NumSol; Type FemEquation;
+ Quantity {
+ { Name u_num; Type Local; NameOfSpace H_num; }
+If(BND_TYPE == BND_Pade)
+If(CRN_TYPE == CRN_Regularization)
+ For m In {1:nPade}
+ { Name u~{m}; Type Local; NameOfSpace H~{m}; }
+ EndFor
+EndIf
+If(CRN_TYPE == CRN_Compatibility)
+ For m In {1:nPade}
+ { Name u~{m}~{0}~{0}; Type Local; NameOfSpace H~{m}~{0}~{0}; }
+ { Name u~{0}~{m}~{0}; Type Local; NameOfSpace H~{0}~{m}~{0}; }
+ { Name u~{0}~{0}~{m}; Type Local; NameOfSpace H~{0}~{0}~{m}; }
+ EndFor
+ For m In {1:nPade}
+ For n In {1:nPade}
+ { Name u~{0}~{m}~{n}; Type Local; NameOfSpace H~{0}~{m}~{n}; }
+ { Name u~{n}~{0}~{m}; Type Local; NameOfSpace H~{n}~{0}~{m}; }
+ { Name u~{m}~{n}~{0}; Type Local; NameOfSpace H~{m}~{n}~{0}; }
+ EndFor
+ EndFor
+ For m In {1:nPade}
+ For n In {1:nPade}
+ For o In {1:nPade}
+ { Name u~{m}~{n}~{o}; Type Local; NameOfSpace H~{m}~{n}~{o}; }
+ EndFor
+ EndFor
+ EndFor
+EndIf
+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; }
+
+// Sommerfeld ABC
+
+If(BND_TYPE == BND_Sommerfeld)
+ Galerkin{ [ -I[]*k[] * Dof{u_num} , {u_num} ]; In FacesAll; Jacobian JSur; Integration I1; }
+EndIf
+
+// Second-order ABC
+
+If(BND_TYPE == BND_Second)
+ Galerkin { [ -I[]*k[] * Dof{u_num} , {u_num} ]; In FacesAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -1/(2*I[]*k[]) * Dof{d u_num} , {d u_num} ]; In FacesAll; Jacobian JSur; Integration I1; }
+If(CRN_TYPE == CRN_Compatibility)
+ Galerkin { [ 3/4 * Dof{u_num} , {u_num} ]; In EdgesAll; Jacobian JLin; Integration I1; }
+ Galerkin { [ 1/(2*I[]*k[]*2*I[]*k[]) * Dof{d u_num} , {d u_num} ]; In EdgesAll; Jacobian JLin; Integration I1; }
+ Galerkin { [ -1/(2*I[]*k[]) * Dof{u_num} , {u_num} ]; In CornersAll; Jacobian JVol; Integration I1; }
+EndIf
+EndIf
+
+// Pade ABC
+
+If(BND_TYPE == BND_Pade)
+If(CRN_TYPE == CRN_Regularization)
+
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In FacesAll; Jacobian JSur; Integration I1; }
+
+ For m In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{m}} , {u_num} ]; In FacesAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u_num} , {u_num} ]; In FacesAll; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u~{m}} , {d u~{m}} ]; In FacesAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{m}+ExpMTheta[]) * Dof{u~{m}} , {u~{m}} ]; In FacesAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{m}+1) * Dof{u_num} , {u~{m}} ]; In FacesAll; Jacobian JSur; Integration I1; }
+ EndFor
+
+EndIf
+If(CRN_TYPE == CRN_Compatibility)
+
+ // --- Faces
+
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In FacesX; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In FacesY; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In FacesZ; Jacobian JSur; Integration I1; }
+
+ For m In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{m}~{0}~{0}} , {u_num} ]; In FacesX; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u_num} , {u_num} ]; In FacesX; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{0}~{m}~{0}} , {u_num} ]; In FacesY; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u_num} , {u_num} ]; In FacesY; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{0}~{0}~{m}} , {u_num} ]; In FacesZ; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u_num} , {u_num} ]; In FacesZ; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u~{m}~{0}~{0}} , {d u~{m}~{0}~{0}} ]; In FacesX; Jacobian JSur; Integration I1; }
+ Galerkin { [ Dof{d u~{0}~{m}~{0}} , {d u~{0}~{m}~{0}} ]; In FacesY; Jacobian JSur; Integration I1; }
+ Galerkin { [ Dof{d u~{0}~{0}~{m}} , {d u~{0}~{0}~{m}} ]; In FacesZ; Jacobian JSur; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+ExpMTheta[]) * Dof{u~{m}~{0}~{0}} , {u~{m}~{0}~{0}} ]; In FacesX; Jacobian JSur; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+ExpMTheta[]) * Dof{u~{0}~{m}~{0}} , {u~{0}~{m}~{0}} ]; In FacesY; Jacobian JSur; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+ExpMTheta[]) * Dof{u~{0}~{0}~{m}} , {u~{0}~{0}~{m}} ]; In FacesZ; Jacobian JSur; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+1) * Dof{u_num} , {u~{m}~{0}~{0}} ]; In FacesX; Jacobian JSur; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+1) * Dof{u_num} , {u~{0}~{m}~{0}} ]; In FacesY; Jacobian JSur; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+1) * Dof{u_num} , {u~{0}~{0}~{m}} ]; In FacesZ; Jacobian JSur; Integration I1; }
+ EndFor
+
+ // --- Edges
+
+ For m In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{0}~{m}~{0}} , {u~{0}~{m}~{0}} ]; In EdgesYZ; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{0}~{0}~{m}} , {u~{0}~{0}~{m}} ]; In EdgesZX; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{m}~{0}~{0}} , {u~{m}~{0}~{0}} ]; In EdgesXY; Jacobian JLin; Integration I1; }
+ EndFor
+ For n In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{0}~{0}~{n}} , {u~{0}~{0}~{n}} ]; In EdgesYZ; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{n}~{0}~{0}} , {u~{n}~{0}~{0}} ]; In EdgesZX; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{0}~{n}~{0}} , {u~{0}~{n}~{0}} ]; In EdgesXY; Jacobian JLin; Integration I1; }
+ EndFor
+
+ For m In {1:nPade}
+ For n In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{0}~{m}~{n}} , {u~{0}~{m}~{0}} ]; In EdgesYZ; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{0}~{m}~{0}} , {u~{0}~{m}~{0}} ]; In EdgesYZ; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{0}~{m}~{n}} , {u~{0}~{0}~{n}} ]; In EdgesYZ; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{0}~{0}~{n}} , {u~{0}~{0}~{n}} ]; In EdgesYZ; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{n}~{0}~{m}} , {u~{0}~{0}~{m}} ]; In EdgesZX; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{0}~{0}~{m}} , {u~{0}~{0}~{m}} ]; In EdgesZX; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{n}~{0}~{m}} , {u~{n}~{0}~{0}} ]; In EdgesZX; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{n}~{0}~{0}} , {u~{n}~{0}~{0}} ]; In EdgesZX; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{m}~{n}~{0}} , {u~{m}~{0}~{0}} ]; In EdgesXY; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{m}~{0}~{0}} , {u~{m}~{0}~{0}} ]; In EdgesXY; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{m}~{n}~{0}} , {u~{0}~{n}~{0}} ]; In EdgesXY; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{0}~{n}~{0}} , {u~{0}~{n}~{0}} ]; In EdgesXY; Jacobian JLin; Integration I1; }
+
+ Galerkin { [ Dof{d u~{0}~{m}~{n}} , {d u~{0}~{m}~{n}} ]; In EdgesYZ; Jacobian JLin; Integration I1; }
+ Galerkin { [ Dof{d u~{n}~{0}~{m}} , {d u~{n}~{0}~{m}} ]; In EdgesZX; Jacobian JLin; Integration I1; }
+ Galerkin { [ Dof{d u~{m}~{n}~{0}} , {d u~{m}~{n}~{0}} ]; In EdgesXY; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+cPade~{n}+ExpMTheta[]) * Dof{u~{0}~{m}~{n}} , {u~{0}~{m}~{n}} ]; In EdgesYZ; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+cPade~{n}+ExpMTheta[]) * Dof{u~{n}~{0}~{m}} , {u~{n}~{0}~{m}} ]; In EdgesZX; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+cPade~{n}+ExpMTheta[]) * Dof{u~{m}~{n}~{0}} , {u~{m}~{n}~{0}} ]; In EdgesXY; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+1) * Dof{u~{0}~{0}~{n}} , {u~{0}~{m}~{n}} ]; In EdgesYZ; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+1) * Dof{u~{n}~{0}~{0}} , {u~{n}~{0}~{m}} ]; In EdgesZX; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+1) * Dof{u~{0}~{n}~{0}} , {u~{m}~{n}~{0}} ]; In EdgesXY; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{n}+1) * Dof{u~{0}~{m}~{0}} , {u~{0}~{m}~{n}} ]; In EdgesYZ; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{n}+1) * Dof{u~{0}~{0}~{m}} , {u~{n}~{0}~{m}} ]; In EdgesZX; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{n}+1) * Dof{u~{m}~{0}~{0}} , {u~{m}~{n}~{0}} ]; In EdgesXY; Jacobian JLin; Integration I1; }
+ EndFor
+ EndFor
+
+ // --- Corners
+
+ For m In {1:nPade}
+ For n In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{0}~{m}~{n}} , {u~{0}~{m}~{n}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{n}~{0}~{m}} , {u~{n}~{0}~{m}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{m}~{n}~{0}} , {u~{m}~{n}~{0}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ EndFor
+ EndFor
+
+ For m In {1:nPade}
+ For n In {1:nPade}
+ For o In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{m}~{n}~{o}} , {u~{0}~{n}~{o}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{0}~{n}~{o}} , {u~{0}~{n}~{o}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{m}~{n}~{o}} , {u~{m}~{0}~{o}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{m}~{0}~{o}} , {u~{m}~{0}~{o}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{o} * Dof{u~{m}~{n}~{o}} , {u~{m}~{n}~{0}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{o} * Dof{u~{m}~{n}~{0}} , {u~{m}~{n}~{0}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ (cPade~{m}+cPade~{n}+cPade~{o}+ExpMTheta[]) * Dof{u~{m}~{n}~{o}} , {u~{m}~{n}~{o}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ (cPade~{m}+1) * Dof{u~{0}~{n}~{o}} , {u~{m}~{n}~{o}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ (cPade~{n}+1) * Dof{u~{m}~{0}~{o}} , {u~{m}~{n}~{o}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ (cPade~{o}+1) * Dof{u~{m}~{n}~{0}} , {u~{m}~{n}~{o}} ]; In CornersAll; Jacobian JSur; Integration I1; }
+ EndFor
+ EndFor
+ EndFor
+
+EndIf
+EndIf
+
+ }
+ }
+}
+
+//==================================================================================================
+// RESOLUTION
+//==================================================================================================
+
+Resolution {
+ { Name NumNormal;
+ System {
+ { Name A; NameOfFormulation NumNormal; Type Real; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A];
+ }
+ }
+ { Name NumCur;
+ System {
+ { Name A; NameOfFormulation NumNormal; Type Real; }
+ { Name B; NameOfFormulation NumCur; Type Real; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A]; PostOperation[NumNormal];
+ Generate[B]; Solve[B]; SaveSolution[B];
+ }
+ }
+ { Name NumSol;
+ System {
+If((BND_TYPE == BND_Pade) && (CRN_TYPE == CRN_Regularization) && (KEPS_TYPE == KEPS_Numeric))
+ { Name A; NameOfFormulation NumNormal; Type Real; }
+ { Name B; NameOfFormulation NumCur; Type Real; }
+EndIf
+ { Name C; NameOfFormulation NumSol; Type Complex; }
+ }
+ Operation {
+If((BND_TYPE == BND_Pade) && (CRN_TYPE == CRN_Regularization) && (KEPS_TYPE == KEPS_Numeric))
+ Generate[A]; Solve[A]; SaveSolution[A]; PostOperation[NumNormal];
+ Generate[B]; Solve[B]; SaveSolution[B]; PostOperation[NumCur];
+EndIf
+ Generate[C]; Solve[C]; SaveSolution[C];
+ }
+ }
+}
+
+//==================================================================================================
+// POSTPRO / POSTOP
+//==================================================================================================
+
+PostProcessing {
+ { Name NumNormal; NameOfFormulation NumNormal;
+ Quantity {
+ { Name u_normal; Value { Local { [ Vector[{u_nx},{u_ny},{u_nz}] / Sqrt[{u_nx}*{u_nx}+{u_ny}*{u_ny}+{u_nz}*{u_nz}] ]; In Region[{FacesAll,BndSca}]; Jacobian JSur; }}}
+ }
+ }
+ { Name NumCur; NameOfFormulation NumCur;
+ Quantity {
+ { Name u_cur; Value { Local { [ {u_cur} ]; In Region[{FacesAll,BndSca}]; Jacobian JSur; }}}
+ }
+ }
+ { Name NumSol; NameOfFormulation NumSol;
+ Quantity {
+ { Name u_ref~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ f_ref[] ]; In Dom; Jacobian JVol; }}}
+ { Name u_num~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ {u_num} ]; In Dom; Jacobian JVol; }}}
+ { Name u_err~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ Norm[f_ref[]-{u_num}]]; In Dom; Jacobian JVol; }}}
+ }
+ }
+ { Name Error; NameOfFormulation NumSol;
+ Quantity {
+ { Name error2; Value { Integral { [ Abs[f_ref[]-{u_num}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[f_ref[]]^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/$energy2Var]] ; In Dom; Jacobian JVol; }}}
+ }
+ }
+}
+
+PostOperation{
+ { Name NumNormal; NameOfPostProcessing NumNormal;
+ Operation {
+ Print [u_normal, OnElementsOf Region[{FacesAll,BndSca}], StoreInField (1001), File "out/u_normal.pos"];
+ }
+ }
+ { Name NumCur; NameOfPostProcessing NumCur;
+ Operation {
+ Print [u_cur, OnElementsOf Region[{FacesAll,BndSca}], StoreInField (1002), File "out/u_cur.pos"];
+ }
+ }
+ { Name NumSol; NameOfPostProcessing NumSol;
+ Operation {
+ tmp1 = Sprintf("out/solRef_cube_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ tmp2 = Sprintf("out/solNum_cube_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ tmp3 = Sprintf("out/solErr_cube_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ Print [u_ref~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Dom, File tmp1];
+ Print [u_num~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Dom, File tmp2];
+ Print [u_err~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Dom, File tmp3];
+ }
+ }
+ { Name Error; NameOfPostProcessing Error;
+ Operation {
+ tmp4 = Sprintf("out/errorAbs_cube_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ tmp5 = Sprintf("out/errorRel_cube_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ Print [error2[Dom], OnRegion Dom, Format Table, StoreInVariable $error2Var];
+ Print [energy2[Dom], OnRegion Dom, Format Table, StoreInVariable $energy2Var];
+ Print [errorAbs, OnRegion Dom, Format Table, SendToServer "Output/1L2-Error (absolute)", File > tmp4];
+ Print [errorRel, OnRegion Dom, Format Table, SendToServer "Output/2L2-Error (relative)", File > tmp5];
+ }
+ }
+}
+
+DefineConstant[
+ R_ = {"NumSol", Name "GetDP/1ResolutionChoices", Visible 1, Choices {"NumNormal", "NumCur", "NumSol"} },
+ P_ = {"NumSol", Name "GetDP/2PostOperationChoices", Visible 1, Choices {"NumNormal", "NumCur", "NumSol", "Error"}}
+];
diff --git a/HelmholtzHABCwithCorners/padeDisk.dat b/HelmholtzHABCwithCorners/padeDisk.dat
new file mode 100644
index 0000000000000000000000000000000000000000..8da1a3ae827d0b3b43dd24fb91e4ddefe4417a95
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padeDisk.dat
@@ -0,0 +1,32 @@
+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;
diff --git a/HelmholtzHABCwithCorners/padeDisk.geo b/HelmholtzHABCwithCorners/padeDisk.geo
new file mode 100644
index 0000000000000000000000000000000000000000..ddc8a090cb58d3af7e9619743d73a16fddf715b0
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padeDisk.geo
@@ -0,0 +1,54 @@
+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};
diff --git a/HelmholtzHABCwithCorners/padeDisk.pro b/HelmholtzHABCwithCorners/padeDisk.pro
new file mode 100644
index 0000000000000000000000000000000000000000..1522df2f5642129eb999b121ca9f02a2ba1e2c4a
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padeDisk.pro
@@ -0,0 +1,314 @@
+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"}}
+];
diff --git a/HelmholtzHABCwithCorners/padePieWedge.dat b/HelmholtzHABCwithCorners/padePieWedge.dat
new file mode 100644
index 0000000000000000000000000000000000000000..940824119befd86a6232b7a563b1efabfa437f51
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padePieWedge.dat
@@ -0,0 +1,46 @@
+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;
diff --git a/HelmholtzHABCwithCorners/padePieWedge.geo b/HelmholtzHABCwithCorners/padePieWedge.geo
new file mode 100644
index 0000000000000000000000000000000000000000..879630efc281be44c8ba659da949e59c998ceed9
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padePieWedge.geo
@@ -0,0 +1,59 @@
+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
+
+
diff --git a/HelmholtzHABCwithCorners/padePieWedge.pro b/HelmholtzHABCwithCorners/padePieWedge.pro
new file mode 100644
index 0000000000000000000000000000000000000000..2bd55e6f0d8551ac96d05c22d96c77724ab1fb5f
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padePieWedge.pro
@@ -0,0 +1,555 @@
+Include "padePieWedge.dat";
+
+//==================================================================================================
+// OPTIONS and PARAMETERS
+//==================================================================================================
+
+BND_Neumann = 0;
+BND_Sommerfeld = 1;
+BND_Second = 2;
+BND_PadeCont = 3;
+BND_PadeDisc = 4;
+CRN_Nothing = 0;
+CRN_Damping = 1;
+CRN_DampingNum = 2;
+CRN_Compatibility = 3;
+CRN_Sommerfeld = 4;
+CRN_DampingNum2 = 5;
+
+DefineConstant[
+ BND_TYPE = {BND_PadeCont,
+ Name "Input/5Model/03Boundary condition (edges)", Highlight "Blue",
+ Choices {BND_Neumann = "Homogeneous Neumann",
+ BND_Sommerfeld = "Sommerfeld ABC",
+ BND_Second = "Second-order ABC",
+ BND_PadeDisc = "Pade ABC (disc at corners)",
+ BND_PadeCont = "Pade ABC (cont at corners)"}},
+ CRN_TYPE = {CRN_Nothing,
+ Name "Input/5Model/04Boundary condition (corners)", Highlight "Red",
+ Visible ((BND_TYPE == BND_Second) || (BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont)),
+ Choices {CRN_Nothing = "Nothing",
+ CRN_Damping = "Damping with keps",
+ CRN_DampingNum = "Damping with keps num",
+ CRN_Compatibility = "Compatibility",
+ CRN_Sommerfeld = "Sommerfeld ABC for Compatibility"}}
+];
+
+DefineConstant[
+ nPade = {4, Min 0, Step 1, Max 6,
+ Name "Input/5Model/05Pade: Number of fields",
+ Visible ((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))},
+ thetaPadeInput = {3, Min 0, Step 1, Max 4,
+ Name "Input/5Model/06Pade: Rotation of branch cut",
+ Visible ((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))}
+];
+
+If((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))
+ If(thetaPadeInput == 0)
+ thetaPade = 0;
+ EndIf
+ If(thetaPadeInput == 1)
+ thetaPade = Pi/8;
+ EndIf
+ If(thetaPadeInput == 2)
+ thetaPade = Pi/4;
+ EndIf
+ If(thetaPadeInput == 3)
+ thetaPade = Pi/3;
+ EndIf
+ If(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
+
+Group {
+ Dom = Region[{DOM}];
+ BndSca = Region[{BND_Scatt}];
+
+If(BND_TYPE == BND_PadeDisc)
+ Edg~{1} = Region[{BND~{1}}];
+ Edg~{2} = Region[{BND~{2}}];
+ Edg~{3} = Region[{BND~{3}}];
+ EdgClo~{1} = Region[{BND~{1},CRN~{1},CRN~{2}}];
+ EdgClo~{2} = Region[{BND~{2},CRN~{2},CRN~{3}}];
+ EdgClo~{3} = Region[{BND~{3},CRN~{3},CRN~{1}}];
+ maxEdg = 3;
+ maxCrn = 3;
+Else
+ Edg~{1} = Region[{BND~{1}}];
+ Edg~{2} = Region[{BND~{2},BND~{3}}];
+ EdgClo~{1} = Region[{BND~{1},CRN~{1},CRN~{2}}];
+ EdgClo~{2} = Region[{BND~{2},BND~{3},CRN~{2},CRN~{1}}];
+ maxEdg = 2;
+ maxCrn = 2;
+EndIf
+ Crn~{1} = Region[{CRN~{1}}];
+ Crn~{2} = Region[{CRN~{2}}];
+ Crn~{3} = Region[{CRN~{3}}];
+
+ CrnAll = Region[{CRN~{1},CRN~{2},CRN~{3}}];
+ EdgAll = Region[{BND~{1},BND~{2},BND~{3}}];
+ DomAll = Region[{Dom,BndSca,EdgAll,CrnAll}];
+
+If(FLAG_REF == REF_Ana)
+ CurvAll = Region[{EdgAll}];
+ DomRef = Region[{DOM}];
+ DomRefAll = Region[{DomRef,BndSca,EdgAll}];
+Else
+ BndExt = Region[{BND_EXT}];
+ BndExtAll = Region[{BND~{1},BND_EXT}];
+ CurvAll = Region[{EdgAll,BndExt}];
+ DomRef = Region[{DOM,DOM_EXT}];
+ DomRefAll = Region[{DomRef,BndSca,BndExtAll}];
+EndIf
+}
+
+Function {
+ RadiusDom[] = XYZ[] + Vector[X_SCA,Y_SCA,0.];
+ NormalGeo[Region[{BND~{1}}]] = RadiusDom[] / Norm[RadiusDom[]];
+ NormalGeo[Region[{BND~{2}}]] = Vector[Cos[alphaDom],-Sin[alphaDom],0.];
+ NormalGeo[Region[{BND~{3}}]] = Vector[-1.,0.,0.];
+If(FLAG_REF == REF_Num)
+ NormalGeo[Region[{BND_EXT}]] = RadiusDom[] / Norm[RadiusDom[]];
+EndIf
+ NormalGeo[BndSca] = XYZ[] / Norm[XYZ[]];
+
+ NumNormal[] = VectorField[XYZ[]]{1001};
+ NumCurv[] = ScalarField[XYZ[]]{1002};
+
+ CurvGeo[Edg~{1}] = 1./dimL;
+ CurvGeo[Edg~{2}] = 0;
+If(BND_TYPE == BND_PadeDisc)
+ CurvGeo[Edg~{3}] = 0;
+EndIf
+If(FLAG_REF == REF_Num)
+ CurvGeo[BndExt] = 1./dimL;
+EndIf
+
+ DistCorner[] = Sqrt[(X[]-X~{3})^2 + (Y[]-Y~{3})^2];
+If(CRN_TYPE == CRN_Damping)
+ //CurvCorner = (2)^(1/2)/LC;
+ //CurvCorner = 2.*Cos[alphaDom/2.]/LC;
+ CurvCorner = 1/(Tan[alphaDom/2.]*LC);
+ Curv[CurvAll] = (DistCorner[] < LC) ? CurvCorner : CurvGeo[] ;
+ CurvEps[CurvAll] = (DistCorner[] < LC) ? CurvCorner : CurvGeo[] ;
+ElseIf(CRN_TYPE == CRN_DampingNum)
+ Curv[CurvAll] = (DistCorner[] < 10*LC) ? NumCurv[] : CurvGeo[] ;
+ CurvEps[CurvAll] = (DistCorner[] < 10*LC) ? NumCurv[] : CurvGeo[] ;
+ElseIf(CRN_TYPE == CRN_DampingNum2)
+ Curv[CurvAll] = CurvGeo[];
+ CurvEps[CurvAll] = (DistCorner[] < 10*LC) ? NumCurv[] : CurvGeo[] ;
+Else
+ Curv[CurvAll] = CurvGeo[];
+ CurvEps[CurvAll] = CurvGeo[];
+EndIf
+
+ kEps[CurvAll] = WAVENUMBER + I[] * 0.39 * WAVENUMBER^(1/3) * (CurvEps[])^(2/3);
+
+ alphaBT[CurvAll] = Curv[]/2 + (Curv[])^2 * 1/(8.*(I[]*k[] - Curv[]));
+ betaBT[CurvAll] = - 1/(2.*(I[]*k[] - Curv[]));
+
+ alphaBTPade[CurvAll] = Curv[]/2 + (Curv[])^2 * 1/(8.*(I[]*k[] - Curv[]));
+ betaBTPade[CurvAll] = - Curv[]/(2*k[]*k[]);
+
+If((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))
+ ExpPTheta[] = Complex[Cos[ thetaPade],Sin[ thetaPade]];
+ ExpMTheta[] = Complex[Cos[-thetaPade],Sin[-thetaPade]];
+ ExpPTheta2[] = Complex[Cos[thetaPade/2.],Sin[thetaPade/2.]];
+EndIf
+}
+
+//==================================================================================================
+// FONCTION SPACES with CONSTRAINTS
+//==================================================================================================
+
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+Constraint {
+ { Name DirichletBC; Case {{ Region BndSca; Value -f_inc[]; }}}
+}
+EndIf
+
+FunctionSpace {
+ { Name H_nx; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{CurvAll}]; Entity NodesOf[All]; }}}
+ { Name H_ny; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{CurvAll}]; Entity NodesOf[All]; }}}
+ { Name H_cur; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{CurvAll}]; Entity NodesOf[All]; }}}
+ { Name H_proj; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{Dom}]; Entity NodesOf[All]; }}}
+If(FLAG_REF == REF_Num)
+ { Name H_ref; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{DomRefAll}]; Entity NodesOf[All]; }}
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+ Constraint {{ NameOfCoef pi; EntityType NodesOf; NameOfConstraint DirichletBC; }}
+EndIf
+ }
+ For i In {1:nPade}
+ { Name H_ref~{i}; Type Form0;BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{BndExtAll}]; Entity NodesOf[All]; }}}
+ EndFor
+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)
+ Constraint {{ NameOfCoef pi; EntityType NodesOf; NameOfConstraint DirichletBC; }}
+EndIf
+ }
+If((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))
+ For iEdg In{1:maxEdg}
+ For i In {1:nPade}
+ { Name H~{iEdg}~{i}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{EdgClo~{iEdg}}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+ For iCrn In{1:maxCrn}
+ For i In {1:nPade}
+ For j In {1:nPade}
+ { Name H~{iCrn}~{i}~{j}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{Crn~{iCrn}}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+ EndFor
+EndIf
+}
+
+//==================================================================================================
+// FORMULATIONS
+//==================================================================================================
+
+Formulation {
+ { Name NumNormal; Type FemEquation;
+ Quantity {
+ { Name u_nx; Type Local; NameOfSpace H_nx; }
+ { Name u_ny; Type Local; NameOfSpace H_ny; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_nx} , {u_nx} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ny} , {u_ny} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompX[NormalGeo[]] , {u_nx} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompY[NormalGeo[]] , {u_ny} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+ }
+ }
+
+ { Name NumCur; Type FemEquation;
+ Quantity {
+ { Name u_nx; Type Local; NameOfSpace H_nx; }
+ { Name u_ny; Type Local; NameOfSpace H_ny; }
+ { Name u_cur; Type Local; NameOfSpace H_cur; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_nx} , {u_nx} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ny} , {u_ny} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompX[NumNormal[]] , {u_nx} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompY[NumNormal[]] , {u_ny} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+
+ Galerkin{ [ Dof{u_cur} , {u_cur} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[1,0,0] * Dof{d u_nx} , {u_cur} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[0,1,0] * Dof{d u_ny} , {u_cur} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+// Galerkin{ [ -CompX[(Vector[0,0,1] /\ NormalGeo[])] * (Vector[0,0,1] /\ NormalGeo[]) * Dof{d u_nx} , {u_cur} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+// Galerkin{ [ -CompY[(Vector[0,0,1] /\ NormalGeo[])] * (Vector[0,0,1] /\ NormalGeo[]) * Dof{d u_ny} , {u_cur} ]; In Region[{CurvAll}]; Jacobian JSur; Integration I1; }
+ }
+ }
+
+If(FLAG_REF == REF_Num)
+ { Name NumRef; Type FemEquation;
+ Quantity {
+ { Name u_ref; Type Local; NameOfSpace H_ref; }
+ For i In {1:nPade}
+ { Name u_ref~{i}; Type Local; NameOfSpace H_ref~{i}; }
+ EndFor
+ }
+ Equation {
+
+ Galerkin{ [ Dof{d u_ref} , {d u_ref} ]; In DomRef; Jacobian JVol; Integration I1; }
+ Galerkin{ [ -k[]^2*Dof{u_ref} , {u_ref} ]; In DomRef; Jacobian JVol; Integration I1; }
+If(FLAG_SIGNAL_BC == SIGNAL_Neumann)
+ Galerkin{ [ -df_inc[] , {u_ref} ]; In BndSca; Jacobian JSur; Integration I1; }
+EndIf
+
+ Galerkin { [ alphaBTPade[] * Dof{u_ref} , {u_ref} ]; In BndExtAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ betaBTPade[] * Dof{d u_ref} , {d u_ref} ]; In BndExtAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_ref} , {u_ref} ]; In BndExtAll; Jacobian JSur; Integration I1; }
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u_ref~{i}} , {u_ref} ]; In BndExtAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u_ref} , {u_ref} ]; In BndExtAll; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u_ref~{i}} , {d u_ref~{i}} ]; In BndExtAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*(ExpPTheta[]*cPade~{i}+1) * Dof{u_ref~{i}} , {u_ref~{i}} ]; In BndExtAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{i}+1) * Dof{u_ref} , {u_ref~{i}} ]; In BndExtAll; Jacobian JSur; Integration I1; }
+ EndFor
+
+ }
+ }
+EndIf
+
+ { Name NumSol; Type FemEquation;
+ Quantity {
+If(FLAG_REF == REF_Num)
+ { Name u_ref; Type Local; NameOfSpace H_ref; }
+EndIf
+ { Name u_num; Type Local; NameOfSpace H_num; }
+If((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))
+ For iEdg In{1:maxEdg}
+ For i In {1:nPade}
+ { Name u~{iEdg}~{i}; Type Local; NameOfSpace H~{iEdg}~{i}; }
+ EndFor
+ EndFor
+ For iCrn In{1:maxCrn}
+ For i In {1:nPade}
+ For j In {1:nPade}
+ { Name u~{iCrn}~{i}~{j}; Type Local; NameOfSpace H~{iCrn}~{i}~{j}; }
+ EndFor
+ EndFor
+ 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
+
+// Sommerfeld ABC
+ If(BND_TYPE == BND_Sommerfeld)
+ Galerkin{ [ -I[]*k[]*Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ EndIf
+
+// Second-order ABC
+ If(BND_TYPE == BND_Second)
+ Galerkin { [ - I[]*k[] * Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ alphaBT[] * Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ betaBT[] * Dof{d u_num} , {d u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ If(CRN_TYPE == CRN_Compatibility)
+ Galerkin { [ 3./4. * Dof{u_num} , {u_num} ]; In CrnAll; Jacobian JLin; Integration I1; }
+ EndIf
+ EndIf
+
+// HABC (auxiliary fields continuous/discontinuous at the corners)
+If((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))
+ For iEdg In{1:maxEdg}
+ Galerkin { [ alphaBTPade[] * Dof{u_num} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ betaBTPade[] * Dof{d u_num} , {d u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u~{iEdg}~{i}} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u_num} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u~{iEdg}~{i}} , {d u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*(ExpPTheta[]*cPade~{i}+1) * Dof{u~{iEdg}~{i}} , {u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{i}+1) * Dof{u_num} , {u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ EndFor
+ EndFor
+
+ For iCrn In{1:maxCrn}
+ iEdg1 = (iCrn == 1) ? maxCrn : iCrn-1;
+ iEdg2 = iCrn;
+ If((iCrn < 3) || (CRN_TYPE == CRN_Compatibility))
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+
+ For j In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iCrn}~{i}~{j}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iCrn}~{j}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+
+ Galerkin { [ (cPade~{i}+cPade~{j}+ExpMTheta[]) * Dof{u~{iCrn}~{i}~{j}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ (cPade~{j}+1) * Dof{u~{iEdg1}~{i}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ (cPade~{i}+1) * Dof{u~{iEdg2}~{j}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ EndFor
+ EndFor
+ EndIf
+ If((iCrn == 3) && (CRN_TYPE == CRN_Sommerfeld))
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[] * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; } // *ExpPTheta2[]
+ Galerkin { [ -I[]*k[] * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; } // *ExpPTheta2[]
+ EndFor
+ EndIf
+ EndFor
+
+EndIf
+ }
+ }
+
+ { Name ProjSol; Type FemEquation;
+ Quantity {
+ { Name u_proj; Type Local; NameOfSpace H_proj; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_proj} , {u_proj} ]; In Dom; Jacobian JVol; Integration I1; }
+ Galerkin{ [ -f_ref[] , {u_proj} ]; In Dom; Jacobian JVol; Integration I1; }
+ }
+ }
+}
+
+//==================================================================================================
+// RESOLUTION
+//==================================================================================================
+
+Resolution {
+ { Name NumNormal;
+ System {
+ { Name A; NameOfFormulation NumNormal; Type Real; NameOfMesh "mainCurv.msh"; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A];
+ }
+ }
+ { Name NumCur;
+ System {
+ { Name A; NameOfFormulation NumNormal; Type Real; NameOfMesh "mainCurv.msh"; }
+ { Name B; NameOfFormulation NumCur; Type Real; NameOfMesh "mainCurv.msh"; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A]; PostOperation[NumNormal];
+ Generate[B]; Solve[B]; SaveSolution[B];
+ }
+ }
+ { Name NumSol;
+ System {
+If((CRN_TYPE == CRN_DampingNum) || (CRN_TYPE == CRN_DampingNum2))
+ { Name A; NameOfFormulation NumNormal; Type Real; NameOfMesh "mainCurv.msh"; }
+ { Name B; NameOfFormulation NumCur; Type Real; NameOfMesh "mainCurv.msh"; }
+EndIf
+If(FLAG_REF == REF_Num)
+ { Name C; NameOfFormulation NumRef; Type Complex; }
+EndIf
+ { Name D; NameOfFormulation NumSol; Type Complex; }
+ }
+ Operation {
+If((CRN_TYPE == CRN_DampingNum) || (CRN_TYPE == CRN_DampingNum2))
+ Generate[A]; Solve[A]; SaveSolution[A]; PostOperation[NumNormal];
+ Generate[B]; Solve[B]; SaveSolution[B]; PostOperation[NumCur];
+EndIf
+If(FLAG_REF == REF_Num)
+ Generate[C]; Solve[C]; SaveSolution[C];
+EndIf
+ Generate[D]; Solve[D]; SaveSolution[D];
+ }
+ }
+ { Name ProjSol;
+ System {{ Name A; NameOfFormulation ProjSol; Type Complex; }}
+ Operation { Generate[A]; Solve[A]; SaveSolution[A]; }
+ }
+}
+
+//==================================================================================================
+// POSTPRO / POSTOP
+//==================================================================================================
+
+PostProcessing {
+ { Name NumNormal; NameOfFormulation NumNormal;
+ Quantity {
+ { Name u_normal; Value { Local { [ Vector[{u_nx},{u_ny},0] / Sqrt[{u_nx}*{u_nx}+{u_ny}*{u_ny}] ]; In Region[{CurvAll}]; Jacobian JSur; }}}
+ }
+ }
+ { Name NumCur; NameOfFormulation NumCur;
+ Quantity {
+ { Name u_cur; Value { Local { [ {u_cur} ]; In Region[{CurvAll}]; Jacobian JSur; }}}
+ }
+ }
+ { Name NumSol; NameOfFormulation NumSol;
+ Quantity {
+If(FLAG_REF == REF_Ana)
+ { Name u_refFull~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}; Value { Local { [ f_ref[] ]; In DomRef; Jacobian JVol; }}}
+ { Name u_ref~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}; Value { Local { [ f_ref[] ]; In Dom; Jacobian JVol; }}}
+ { Name u_num~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}; Value { Local { [ {u_num} ]; In Dom; Jacobian JVol; }}}
+ { Name u_err~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}; Value { Local { [ f_ref[]-{u_num} ]; In Dom; Jacobian JVol; }}}
+ElseIf(FLAG_REF == REF_Num)
+ { Name u_refFull~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}; Value { Local { [ {u_ref} ]; In DomRef; Jacobian JVol; }}}
+ { Name u_ref~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}; Value { Local { [ {u_ref} ]; In Dom; Jacobian JVol; }}}
+ { Name u_num~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}; Value { Local { [ {u_num} ]; In Dom; Jacobian JVol; }}}
+ { Name u_err~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}; Value { Local { [ {u_ref}-{u_num} ]; In Dom; Jacobian JVol; }}}
+EndIf
+ For iEdg In{1:maxEdg}
+ For i In {1:nPade}
+ { Name u_numAux~{iEdg}~{i}~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}; Value { Local { [ {u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; }}}
+ EndFor
+ EndFor
+ }
+ }
+ { Name Errors; NameOfFormulation NumSol;
+ Quantity {
+If(FLAG_REF == REF_Ana)
+ { Name error2; Value { Integral { [ Norm[f_ref[]-{u_num}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[f_ref[]]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ElseIf(FLAG_REF == REF_Num)
+ { Name error2; Value { Integral { [ Norm[{u_ref}-{u_num}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[{u_ref}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+EndIf
+ { Name errorAbs; Value { Term { Type Global; [Sqrt[$error2Var]] ; In Dom; Jacobian JVol; }}}
+ { Name errorRel; Value { Term { Type Global; [Sqrt[$error2Var/$energy2Var]] ; In Dom; Jacobian JVol; }}}
+ }
+ }
+ { Name ProjErrors; NameOfFormulation ProjSol;
+ Quantity {
+ { Name error2; Value { Integral { [ Norm[f_ref[]-{u_proj}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[f_ref[]]^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/$energy2Var]] ; In Dom; Jacobian JVol; }}}
+ }
+ }
+}
+
+PostOperation{
+ { Name NumNormal; NameOfPostProcessing NumNormal;
+ Operation {
+ tmp1 = Sprintf("out/u_normal_%g_%g.pos", alphaDomSave, LC);
+ Print [u_normal, OnElementsOf Region[{BND~{2}}], File tmp1];
+ Print [u_normal, OnElementsOf Region[{CurvAll}], StoreInField (1001)];
+ }
+ }
+ { Name NumCur; NameOfPostProcessing NumCur;
+ Operation {
+ tmp1 = Sprintf("out/u_cur_%g_%g.pos", alphaDomSave, LC);
+ Print [u_cur, OnElementsOf Region[{BND~{2}}], File tmp1];
+ Print [u_cur, OnElementsOf Region[{CurvAll}], StoreInField (1002)];
+ }
+ }
+ { Name NumSol; NameOfPostProcessing NumSol;
+ Operation {
+ tmp1 = Sprintf("out/solRefFull_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, alphaDomSave);
+ tmp2 = Sprintf("out/solRef_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, alphaDomSave);
+ tmp3 = Sprintf("out/solNum_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, alphaDomSave);
+ tmp4 = Sprintf("out/solErr_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, alphaDomSave);
+// Print [u_refFull~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}, OnElementsOf DomRef, File tmp1];
+ Print [u_ref~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}, OnElementsOf Dom, File tmp2];
+ Print [u_num~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}, OnElementsOf Dom, File tmp3];
+ Print [u_err~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}, OnElementsOf Dom, File tmp4];
+// For iEdg In{1:maxEdg}
+// For i In {1:nPade}
+// tmp5 = Sprintf("out/solNum_%g_%g_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, alphaDomSave, iEdg, i);
+// Print [u_numAux~{iEdg}~{i}~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{alphaDomSave}, OnElementsOf Edg~{iEdg}, File tmp5];
+// EndFor
+// EndFor
+ }
+ }
+ { Name Errors; NameOfPostProcessing Errors;
+ Operation {
+ tmp5 = Sprintf("out/errorAbs_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ tmp6 = Sprintf("out/errorRel_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ Print [error2[Dom], OnRegion Dom, Format Table, StoreInVariable $error2Var];
+ Print [energy2[Dom], OnRegion Dom, Format Table, StoreInVariable $energy2Var];
+ Print [errorAbs, OnRegion Dom, Format Table, SendToServer "Output/1L2-Error (absolute)", File > tmp5];
+ Print [errorRel, OnRegion Dom, Format Table, SendToServer "Output/2L2-Error (relative)", File > tmp6];
+ }
+ }
+ { Name ProjErrors; NameOfPostProcessing ProjErrors;
+ Operation {
+ tmp1 = Sprintf("out/errorAbs_Proj.dat");
+ tmp2 = Sprintf("out/errorRel_Proj.dat");
+ Print [error2[Dom], OnRegion Dom, Format Table, StoreInVariable $error2Var];
+ Print [energy2[Dom], OnRegion Dom, Format Table, StoreInVariable $energy2Var];
+ Print [errorAbs, OnRegion Dom, Format Table, SendToServer "Output/1L2-Error (absolute)", File > tmp1];
+ Print [errorRel, OnRegion Dom, Format Table, SendToServer "Output/2L2-Error (relative)", File > tmp2];
+ }
+ }
+}
+
+DefineConstant[
+ R_ = {"NumSol", Name "GetDP/1ResolutionChoices", Visible 1, Choices {"NumNormal", "NumCur", "NumSol", "ProjSol"} },
+ P_ = {"NumSol, Errors", Name "GetDP/2PostOperationChoices", Visible 1, Choices {"NumNormal", "NumCur", "NumSol", "Errors", "ProjErrors"}}
+];
diff --git a/HelmholtzHABCwithCorners/padePolygon.dat b/HelmholtzHABCwithCorners/padePolygon.dat
new file mode 100644
index 0000000000000000000000000000000000000000..719d4bd84937c1043411de3a7a9a0586038ae05e
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padePolygon.dat
@@ -0,0 +1,37 @@
+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;
diff --git a/HelmholtzHABCwithCorners/padePolygon.geo b/HelmholtzHABCwithCorners/padePolygon.geo
new file mode 100644
index 0000000000000000000000000000000000000000..a78f4d59d1ac2ab365bd72e868081eb828a7ff46
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padePolygon.geo
@@ -0,0 +1,77 @@
+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
diff --git a/HelmholtzHABCwithCorners/padePolygon.pro b/HelmholtzHABCwithCorners/padePolygon.pro
new file mode 100644
index 0000000000000000000000000000000000000000..01f7839907d6b5eb76b580702953163e4334482f
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padePolygon.pro
@@ -0,0 +1,630 @@
+Include "padePolygon.dat";
+
+//==================================================================================================
+// OPTIONS and PARAMETERS
+//==================================================================================================
+
+BND_Neumann = 0;
+BND_Sommerfeld = 1;
+BND_Second = 2;
+BND_PadeCont = 3;
+BND_PadeDisc = 4;
+CRN_Nothing = 0;
+CRN_Damping = 1;
+CRN_DampingNum = 2;
+CRN_Compatibility = 3;
+CRN_Sommerfeld = 4;
+CRN_DampingNum2 = 5;
+
+DefineConstant[
+ BND_TYPE = {BND_PadeCont,
+ Name "Input/5Model/03Boundary condition (edges)", Highlight "Blue",
+ Choices {BND_Neumann = "Homogeneous Neumann",
+ BND_Sommerfeld = "Sommerfeld ABC",
+ BND_Second = "Second-order ABC",
+ BND_PadeDisc = "Pade ABC (disc at corners)",
+ BND_PadeCont = "Pade ABC (cont at corners)"}},
+ CRN_TYPE = {CRN_DampingNum,
+ Name "Input/5Model/04Boundary condition (corners)", Highlight "Red",
+ Visible ((BND_TYPE == BND_Second) || (BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont)),
+ Choices {CRN_Nothing = "Nothing",
+ CRN_DampingNum2 = "Damping with num curv (without aux terms)",
+ CRN_DampingNum = "Damping with num curv",
+ CRN_Damping = "Damping with formula curv",
+ CRN_Compatibility = "Compatibility",
+ CRN_Sommerfeld = "Sommerfeld ABC at Corners"}}
+];
+
+DefineConstant[
+ nPade = {4, Choices {0, 1, 2, 3, 4, 5, 6, 7, 8},
+ Name "Input/5Model/05Pade: Number of fields",
+ Visible ((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))},
+ thetaPadeInput = {3, Min 0, Step 1, Max 4,
+ Name "Input/5Model/06Pade: Rotation of branch cut",
+ Visible ((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))}
+];
+
+If((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))
+ If(thetaPadeInput == 0)
+ thetaPade = 0;
+ EndIf
+ If(thetaPadeInput == 1)
+ thetaPade = Pi/8;
+ EndIf
+ If(thetaPadeInput == 2)
+ thetaPade = Pi/4;
+ EndIf
+ If(thetaPadeInput == 3)
+ thetaPade = Pi/3;
+ EndIf
+ If(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
+
+Group {
+ Dom = Region[{DOM}];
+ BndSca = Region[{BND_SCA}];
+
+ CrnAll = {};
+ EdgAll = {};
+ For iCrn In{1:DOM_Num}
+ Crn~{iCrn} = Region[{CRN~{iCrn}}];
+ CrnAll += Region[{CRN~{iCrn}}];
+ EndFor
+ For iEdg In{1:DOM_Num}
+ iCrn1 = (iEdg == 1) ? DOM_Num : iEdg-1;
+ iCrn2 = iEdg;
+ Edg~{iEdg} = Region[{BND~{iEdg}}];
+ EdgClo~{iEdg} = Region[{BND~{iEdg},CRN~{iCrn1},CRN~{iCrn2}}];
+ EdgAll += Region[{BND~{iEdg}}];
+ EndFor
+
+ DomAll = Region[{Dom,BndSca,EdgAll,CrnAll}];
+
+If(FLAG_REF == REF_Ana)
+ DomRef = Region[{DOM}];
+ DomRefAll = Region[{DomRef,BndSca,EdgAll,CrnAll}];
+Else
+ CrnRefAll = {};
+ EdgRefAll = {};
+ For iCrn In{1:4}
+ CrnRef~{iCrn} = Region[{CRN_EXT~{iCrn}}];
+ CrnRefAll += Region[{CRN_EXT~{iCrn}}];
+ EndFor
+ For iEdg In{1:4}
+ iCrn1 = (iEdg == 1) ? 4 : iEdg-1;
+ iCrn2 = iEdg;
+ EdgRef~{iEdg} = Region[{BND_EXT~{iEdg}}];
+ EdgRefClo~{iEdg} = Region[{BND_EXT~{iEdg},CRN_EXT~{iCrn1},CRN_EXT~{iCrn2}}];
+ EdgRefAll += Region[{BND_EXT~{iEdg}}];
+ EndFor
+
+ DomRef = Region[{DOM,DOM_EXT}];
+ DomRefAll = Region[{DomRef,BndSca,EdgRefAll,CrnRefAll}];
+EndIf
+}
+
+Function {
+
+ NormalNum[] = VectorField[XYZ[]]{1001};
+ CurvNum[] = ScalarField[XYZ[]]{1002};
+
+ For iEdg In{1:DOM_Num}
+ AngleNormal = (iEdg-1)*angleInterior*Pi/180 - Pi/2;
+ NormalGeo[Edg~{iEdg}] = Vector[Cos[AngleNormal], Sin[AngleNormal], 0.];
+ EndFor
+
+If(CRN_TYPE == CRN_Damping)
+ //CurvCorner = (2)^(1/2)/LC;
+ alphaDom = Pi - 2*Pi/DOM_Num;
+ CurvCorner = 1/(Tan[alphaDom/2.]*LC);
+ distMin~{0}[] = L_EXT;
+For i In{1:DOM_Num}
+ dist~{i}[] = Sqrt[(X[]-X~{i})^2 + (Y[]-Y~{i})^2];
+ distMin~{i}[] = (distMin~{i-1}[] < dist~{i}[]) ? distMin~{i-1}[] : dist~{i}[];
+EndFor
+ Curv[EdgAll] = (distMin~{DOM_Num}[] < LC) ? CurvCorner : 0;
+ElseIf((CRN_TYPE == CRN_DampingNum) || (CRN_TYPE == CRN_DampingNum2))
+If(DOM_Num == 1)
+ Curv[EdgAll] = 1/DOM_R;
+Else
+ Curv[EdgAll] = CurvNum[];
+EndIf
+Else
+ Curv[EdgAll] = 0.;
+EndIf
+
+ kEps[EdgAll] = WAVENUMBER + I[] * 0.39 * WAVENUMBER^(1/3) * (Curv[])^(2/3);
+ alphaBTPade[EdgAll] = Curv[]/2 + (Curv[])^2 * 1/(8.*(I[]*k[] - Curv[]));
+ betaBTPade[EdgAll] = - Curv[]/(2*k[]*k[]);
+
+If((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))
+ ExpPTheta[] = Complex[Cos[ thetaPade],Sin[ thetaPade]];
+ ExpMTheta[] = Complex[Cos[-thetaPade],Sin[-thetaPade]];
+ ExpPTheta2[] = Complex[Cos[thetaPade/2.],Sin[thetaPade/2.]];
+EndIf
+
+If(FLAG_REF == REF_Num)
+ nPadeRef = 8;
+ mPadeRef = 2*nPadeRef+1;
+ thetaPadeRef = Pi/3;
+ For j In{1:nPadeRef}
+ cPadeRef~{j} = Tan[j*Pi/mPadeRef]^2;
+ EndFor
+ ExpPThetaRef[] = Complex[Cos[ thetaPadeRef],Sin[ thetaPadeRef]];
+ ExpMThetaRef[] = Complex[Cos[-thetaPadeRef],Sin[-thetaPadeRef]];
+ ExpPThetaRef2[] = Complex[Cos[thetaPadeRef/2.],Sin[thetaPadeRef/2.]];
+EndIf
+}
+
+//==================================================================================================
+// FONCTION SPACES with CONSTRAINTS
+//==================================================================================================
+
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+Constraint {
+ { Name DirichletBC; Case {{ Region BndSca; Value -f_inc[]; }}}
+}
+EndIf
+
+FunctionSpace {
+ { Name H_nx; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{EdgAll}]; Entity NodesOf[All]; }}}
+ { Name H_ny; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{EdgAll}]; Entity NodesOf[All]; }}}
+ { Name H_cur; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{EdgAll}]; Entity NodesOf[All]; }}}
+ { Name H_proj; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{Dom}]; Entity NodesOf[All]; }}}
+If(FLAG_REF == REF_Num)
+ { Name H_ref; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{DomRefAll}]; Entity NodesOf[All]; }}
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+ Constraint {{ NameOfCoef pi; EntityType NodesOf; NameOfConstraint DirichletBC; }}
+EndIf
+ }
+ For iEdg In {1:4}
+ For i In {1:nPadeRef}
+ { Name H_ref~{iEdg}~{i}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{EdgRefClo~{iEdg}}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+ For iCrn In {1:4}
+ For i In {1:nPadeRef}
+ For j In {1:nPadeRef}
+ { Name H_ref~{iCrn}~{i}~{j}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{CrnRef~{iCrn}}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+ EndFor
+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)
+ Constraint {{ NameOfCoef pi; EntityType NodesOf; NameOfConstraint DirichletBC; }}
+EndIf
+ }
+If(BND_TYPE == BND_PadeCont)
+ For i In {1:nPade}
+ { Name H~{i}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{EdgAll}]; Entity NodesOf[All]; }}}
+ EndFor
+EndIf
+If(BND_TYPE == BND_PadeDisc)
+ For iEdg In {1:DOM_Num}
+ For i In {1:nPade}
+ { Name H~{iEdg}~{i}; Type Form0;
+ BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{EdgClo~{iEdg}}]; Entity NodesOf[All]; }}
+ }
+ EndFor
+ EndFor
+If(CRN_TYPE == CRN_Compatibility)
+ For iCrn In {1:DOM_Num}
+ For i In {1:nPade}
+ For j In {1:nPade}
+ { Name H~{iCrn}~{i}~{j}; Type Form0;
+ BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{Crn~{iCrn}}]; Entity NodesOf[All]; }}
+ }
+ EndFor
+ EndFor
+ EndFor
+EndIf
+EndIf
+}
+
+//==================================================================================================
+// FORMULATIONS
+//==================================================================================================
+
+Formulation {
+ { Name NumNormal; Type FemEquation;
+ Quantity {
+ { Name u_nx; Type Local; NameOfSpace H_nx; }
+ { Name u_ny; Type Local; NameOfSpace H_ny; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_nx} , {u_nx} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ny} , {u_ny} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompX[NormalGeo[]] , {u_nx} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompY[NormalGeo[]] , {u_ny} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+ }
+ }
+
+ { Name NumCur; Type FemEquation;
+ Quantity {
+ { Name u_nx; Type Local; NameOfSpace H_nx; }
+ { Name u_ny; Type Local; NameOfSpace H_ny; }
+ { Name u_cur; Type Local; NameOfSpace H_cur; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_nx} , {u_nx} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ny} , {u_ny} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompX[NormalNum[]] , {u_nx} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompY[NormalNum[]] , {u_ny} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+
+ Galerkin{ [ Dof{u_cur} , {u_cur} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[1,0,0] * Dof{d u_nx} , {u_cur} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[0,1,0] * Dof{d u_ny} , {u_cur} ]; In Region[{EdgAll}]; Jacobian JSur; Integration I1; }
+ }
+ }
+
+If(FLAG_REF == REF_Num)
+ { Name NumRef; Type FemEquation;
+ Quantity {
+ { Name u_ref; Type Local; NameOfSpace H_ref; }
+ For iEdg In {1:4}
+ For i In {1:nPadeRef}
+ { Name u_ref~{iEdg}~{i}; Type Local; NameOfSpace H_ref~{iEdg}~{i}; }
+ EndFor
+ EndFor
+ For iCrn In {1:4}
+ For i In {1:nPadeRef}
+ For j In {1:nPadeRef}
+ { Name u_ref~{iCrn}~{i}~{j}; Type Local; NameOfSpace H_ref~{iCrn}~{i}~{j}; }
+ EndFor
+ EndFor
+ EndFor
+ }
+ Equation {
+
+ Galerkin{ [ Dof{d u_ref} , {d u_ref} ]; In DomRef; Jacobian JVol; Integration I1; }
+ Galerkin{ [ -k[]^2*Dof{u_ref} , {u_ref} ]; In DomRef; Jacobian JVol; Integration I1; }
+If(FLAG_SIGNAL_BC == SIGNAL_Neumann)
+ Galerkin{ [ -df_inc[] , {u_ref} ]; In BndSca; Jacobian JSur; Integration I1; }
+EndIf
+
+ For iEdg In {1:4}
+ Galerkin { [ -I[]*k[]*ExpPThetaRef2[] * Dof{u_ref} , {u_ref} ]; In EdgRef~{iEdg}; Jacobian JSur; Integration I1; }
+ For i In{1:nPadeRef}
+ Galerkin { [ -I[]*k[]*ExpPThetaRef2[] * 2./mPadeRef * cPadeRef~{i} * Dof{u_ref~{iEdg}~{i}} , {u_ref} ]; In EdgRef~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPThetaRef2[] * 2./mPadeRef * cPadeRef~{i} * Dof{u_ref} , {u_ref} ]; In EdgRef~{iEdg}; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u_ref~{iEdg}~{i}} , {d u_ref~{iEdg}~{i}} ]; In EdgRef~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -k[]^2*(ExpPThetaRef[]*cPadeRef~{i}+1) * Dof{u_ref~{iEdg}~{i}} , {u_ref~{iEdg}~{i}} ]; In EdgRef~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPThetaRef[]*(cPadeRef~{i}+1) * Dof{u_ref} , {u_ref~{iEdg}~{i}} ]; In EdgRef~{iEdg}; Jacobian JSur; Integration I1; }
+ EndFor
+ EndFor
+ For iCrn In {1:4}
+ iEdg1 = iCrn;
+ iEdg2 = (iCrn == 4) ? 1 : iCrn+1;
+ For i In{1:nPadeRef}
+ Galerkin { [ -I[]*k[]*ExpPThetaRef2[] * Dof{u_ref~{iEdg1}~{i}} , {u_ref~{iEdg1}~{i}} ]; In CrnRef~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPThetaRef2[] * Dof{u_ref~{iEdg2}~{i}} , {u_ref~{iEdg2}~{i}} ]; In CrnRef~{iCrn}; Jacobian JLin; Integration I1; }
+ For j In{1:nPadeRef}
+ Galerkin { [ -I[]*k[]*ExpPThetaRef2[] * 2./mPadeRef * cPadeRef~{j} * Dof{u_ref~{iEdg1}~{i}} , {u_ref~{iEdg1}~{i}} ]; In CrnRef~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPThetaRef2[] * 2./mPadeRef * cPadeRef~{j} * Dof{u_ref~{iEdg2}~{i}} , {u_ref~{iEdg2}~{i}} ]; In CrnRef~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPThetaRef2[] * 2./mPadeRef * cPadeRef~{j} * Dof{u_ref~{iCrn}~{i}~{j}} , {u_ref~{iEdg1}~{i}} ]; In CrnRef~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPThetaRef2[] * 2./mPadeRef * cPadeRef~{j} * Dof{u_ref~{iCrn}~{j}~{i}} , {u_ref~{iEdg2}~{i}} ]; In CrnRef~{iCrn}; Jacobian JLin; Integration I1; }
+
+ Galerkin { [ (cPadeRef~{i}+cPadeRef~{j}+ExpMThetaRef[]) * Dof{u_ref~{iCrn}~{i}~{j}} , {u_ref~{iCrn}~{i}~{j}} ]; In CrnRef~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ (cPadeRef~{j}+1) * Dof{u_ref~{iEdg1}~{i}} , {u_ref~{iCrn}~{i}~{j}} ]; In CrnRef~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ (cPadeRef~{i}+1) * Dof{u_ref~{iEdg2}~{j}} , {u_ref~{iCrn}~{i}~{j}} ]; In CrnRef~{iCrn}; Jacobian JLin; Integration I1; }
+ EndFor
+ EndFor
+ EndFor
+
+ }
+ }
+EndIf
+
+ { Name NumSol; Type FemEquation;
+ Quantity {
+If(FLAG_REF == REF_Num)
+ { Name u_ref; Type Local; NameOfSpace H_ref; }
+EndIf
+ { Name u_num; Type Local; NameOfSpace H_num; }
+If(BND_TYPE == BND_PadeCont)
+ For i In {1:nPade}
+ { Name u~{i}; Type Local; NameOfSpace H~{i}; }
+ EndFor
+EndIf
+If(BND_TYPE == BND_PadeDisc)
+ For iEdg In {1:DOM_Num}
+ For i In {1:nPade}
+ { Name u~{iEdg}~{i}; Type Local; NameOfSpace H~{iEdg}~{i}; }
+ EndFor
+ EndFor
+If(CRN_TYPE == CRN_Compatibility)
+ For iCrn In {1:DOM_Num}
+ For i In {1:nPade}
+ For j In {1:nPade}
+ { Name u~{iCrn}~{i}~{j}; Type Local; NameOfSpace H~{iCrn}~{i}~{j}; }
+ EndFor
+ EndFor
+ EndFor
+EndIf
+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
+
+// Sommerfeld ABC
+If(BND_TYPE == BND_Sommerfeld)
+ Galerkin{ [ -I[]*k[]*Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+EndIf
+
+// Second-order ABC
+If(BND_TYPE == BND_Second)
+ Galerkin { [ - I[]*k[] * Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ - 1/(2*I[]*kEps[]) * Dof{d u_num} , {d u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+If(CRN_TYPE == CRN_Compatibility)
+ Galerkin { [ 3./4. * Dof{u_num} , {u_num} ]; In CrnAll; Jacobian JLin; Integration I1; }
+EndIf
+EndIf
+
+// HABC (continuity at corners)
+If(BND_TYPE == BND_PadeCont)
+If((CRN_TYPE == CRN_Damping) || (CRN_TYPE == CRN_DampingNum))
+ Galerkin { [ alphaBTPade[] * Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ betaBTPade[] * Dof{d u_num} , {d u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+EndIf
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u~{i}} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u~{i}} , {d u~{i}} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*(ExpPTheta[]*cPade~{i}+1) * Dof{u~{i}} , {u~{i}} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{i}+1) * Dof{u_num} , {u~{i}} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ EndFor
+EndIf
+
+// HABC (discontinuity at corners)
+
+If(BND_TYPE == BND_PadeDisc)
+ For iEdg In {1:DOM_Num}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u~{iEdg}~{i}} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u_num} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u~{iEdg}~{i}} , {d u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*(ExpPTheta[]*cPade~{i}+1) * Dof{u~{iEdg}~{i}} , {u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{i}+1) * Dof{u_num} , {u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ EndFor
+ EndFor
+
+If(CRN_TYPE == CRN_Sommerfeld)
+ For iCrn In {1:DOM_Num}
+ iEdg1 = iCrn;
+ iEdg2 = (iCrn == DOM_Num) ? 1 : iCrn+1;
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[] * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[] * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ EndFor
+ EndFor
+EndIf
+
+If(CRN_TYPE == CRN_Compatibility)
+ For iCrn In {1:DOM_Num}
+ iEdg1 = iCrn;
+ iEdg2 = (iCrn == DOM_Num) ? 1 : iCrn+1;
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ For j In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iCrn}~{i}~{j}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iCrn}~{j}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+
+ Galerkin { [ (cPade~{i}+cPade~{j}+ExpMTheta[]) * Dof{u~{iCrn}~{i}~{j}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ (cPade~{j}+1) * Dof{u~{iEdg1}~{i}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ (cPade~{i}+1) * Dof{u~{iEdg2}~{j}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ EndFor
+ EndFor
+ EndFor
+EndIf
+EndIf
+
+ }
+ }
+
+ { Name ProjSol; Type FemEquation;
+ Quantity {
+ { Name u_proj; Type Local; NameOfSpace H_proj; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_proj} , {u_proj} ]; In Dom; Jacobian JVol; Integration I1; }
+ Galerkin{ [ -f_ref[] , {u_proj} ]; In Dom; Jacobian JVol; Integration I1; }
+ }
+ }
+}
+
+//==================================================================================================
+// RESOLUTION
+//==================================================================================================
+
+Resolution {
+ { Name NumNormal;
+ System {
+ { Name A; NameOfFormulation NumNormal; Type Real; NameOfMesh "mainCurv.msh"; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A];
+ }
+ }
+ { Name NumCur;
+ System {
+ { Name A; NameOfFormulation NumNormal; Type Real; NameOfMesh "mainCurv.msh"; }
+ { Name B; NameOfFormulation NumCur; Type Real; NameOfMesh "mainCurv.msh"; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A]; PostOperation[NumNormal];
+ Generate[B]; Solve[B]; SaveSolution[B];
+ }
+ }
+ { Name NumSol;
+ System {
+If((CRN_TYPE == CRN_DampingNum) || (CRN_TYPE == CRN_DampingNum2))
+ { Name A; NameOfFormulation NumNormal; Type Real; NameOfMesh "mainCurv.msh"; }
+ { Name B; NameOfFormulation NumCur; Type Real; NameOfMesh "mainCurv.msh"; }
+EndIf
+If(FLAG_REF == REF_Num)
+ { Name C; NameOfFormulation NumRef; Type Complex; }
+EndIf
+ { Name D; NameOfFormulation NumSol; Type Complex; }
+ }
+ Operation {
+If((CRN_TYPE == CRN_DampingNum) || (CRN_TYPE == CRN_DampingNum2))
+ Generate[A]; Solve[A]; SaveSolution[A]; PostOperation[NumNormal];
+ Generate[B]; Solve[B]; SaveSolution[B]; PostOperation[NumCur];
+EndIf
+If(FLAG_REF == REF_Num)
+ Generate[C]; Solve[C]; SaveSolution[C];
+EndIf
+ Generate[D]; Solve[D]; SaveSolution[D];
+ }
+ }
+ { Name ProjSol;
+ System {{ Name A; NameOfFormulation ProjSol; Type Complex; }}
+ Operation { Generate[A]; Solve[A]; SaveSolution[A]; }
+ }
+}
+
+//==================================================================================================
+// POSTPRO / POSTOP
+//==================================================================================================
+
+PostProcessing {
+ { Name NumNormal; NameOfFormulation NumNormal;
+ Quantity {
+ { Name u_normal; Value { Local { [ Vector[{u_nx},{u_ny},0] / Sqrt[{u_nx}*{u_nx}+{u_ny}*{u_ny}] ]; In Region[{EdgAll}]; Jacobian JSur; }}}
+ }
+ }
+ { Name NumCur; NameOfFormulation NumCur;
+ Quantity {
+ { Name u_cur; Value { Local { [ {u_cur} ]; In Region[{EdgAll}]; Jacobian JSur; }}}
+ }
+ }
+ { Name NumSol; NameOfFormulation NumSol;
+ Quantity {
+If(FLAG_REF == REF_Ana)
+ { Name u_ref1~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}; Value { Local { [ f_ref[] ]; In DomRef; Jacobian JVol; }}}
+ { Name u_ref2~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}; Value { Local { [ f_ref[] ]; In Dom; Jacobian JVol; }}}
+ { Name u_num~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}; Value { Local { [ {u_num} ]; In Dom; Jacobian JVol; }}}
+ { Name u_err~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}; Value { Local { [ f_ref[]-{u_num} ]; In Dom; Jacobian JVol; }}}
+ElseIf(FLAG_REF == REF_Num)
+ { Name u_ref1~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}; Value { Local { [ {u_ref} ]; In DomRef; Jacobian JVol; }}}
+ { Name u_ref2~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}; Value { Local { [ {u_ref} ]; In Dom; Jacobian JVol; }}}
+ { Name u_num~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}; Value { Local { [ {u_num} ]; In Dom; Jacobian JVol; }}}
+ { Name u_err~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}; Value { Local { [ {u_ref}-{u_num} ]; In Dom; Jacobian JVol; }}}
+EndIf
+If(FLAG_REF == REF_Ana)
+ { Name error2; Value { Integral { [ Norm[f_ref[]-{u_num}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[f_ref[]]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ElseIf(FLAG_REF == REF_Num)
+ { Name error2; Value { Integral { [ Norm[{u_ref}-{u_num}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[{u_ref}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+EndIf
+ { Name errorAbs; Value { Term { Type Global; [Sqrt[$error2Var]] ; In Dom; Jacobian JVol; }}}
+ { Name errorRel; Value { Term { Type Global; [Sqrt[$error2Var/$energy2Var]] ; In Dom; Jacobian JVol; }}}
+ }
+ }
+ { Name Errors; NameOfFormulation NumSol;
+ Quantity {
+If(FLAG_REF == REF_Ana)
+ { Name error2; Value { Integral { [ Norm[f_ref[]-{u_num}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[f_ref[]]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ElseIf(FLAG_REF == REF_Num)
+ { Name error2; Value { Integral { [ Norm[{u_ref}-{u_num}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[{u_ref}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+EndIf
+ { Name errorAbs; Value { Term { Type Global; [Sqrt[$error2Var]] ; In Dom; Jacobian JVol; }}}
+ { Name errorRel; Value { Term { Type Global; [Sqrt[$error2Var/$energy2Var]] ; In Dom; Jacobian JVol; }}}
+ }
+ }
+ { Name ProjErrors; NameOfFormulation ProjSol;
+ Quantity {
+ { Name u_ref ; Value { Local { [ f_ref[] ]; In Dom; Jacobian JVol; }}}
+ { Name u_proj; Value { Local { [ {u_proj} ]; In Dom; Jacobian JVol; }}}
+ { Name u_diff; Value { Local { [ {u_proj}-f_ref[] ]; In Dom; Jacobian JVol; }}}
+ { Name error2; Value { Integral { [ Norm[f_ref[]-{u_proj}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[f_ref[]]^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/$energy2Var]] ; In Dom; Jacobian JVol; }}}
+ }
+ }
+}
+
+PostOperation{
+ { Name NumNormal; NameOfPostProcessing NumNormal;
+ Operation {
+ Print [u_normal, OnElementsOf Region[{EdgAll}], StoreInField (1001)];
+ }
+ }
+ { Name NumCur; NameOfPostProcessing NumCur;
+ Operation {
+ tmp0 = Sprintf("out/numCur_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ Print [u_cur, OnElementsOf Region[{EdgAll}], StoreInField (1002), File tmp0];
+ }
+ }
+ { Name NumSol; NameOfPostProcessing NumSol;
+ Operation {
+ tmp1 = Sprintf("out/solRef1_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, DOM_Num);
+ tmp2 = Sprintf("out/solRef2_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, DOM_Num);
+ tmp3 = Sprintf("out/solNum_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, DOM_Num);
+ tmp4 = Sprintf("out/solErr_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, DOM_Num);
+ Print [u_ref1~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}, OnElementsOf DomRef, File tmp1];
+ Print [u_ref2~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}, OnElementsOf Dom, File tmp2];
+ Print [u_num~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}, OnElementsOf Dom, File tmp3];
+ Print [u_err~{FLAG_SIGNAL_BC}~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}~{DOM_Num}, OnElementsOf Dom, File tmp4];
+ }
+ }
+ { Name Errors; NameOfPostProcessing Errors;
+ Operation {
+ tmp5 = Sprintf("out/errorAbs_poly_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ tmp6 = Sprintf("out/errorRel_poly_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ Print [error2[Dom], OnRegion Dom, Format Table, StoreInVariable $error2Var];
+ Print [energy2[Dom], OnRegion Dom, Format Table, StoreInVariable $energy2Var];
+ Print [errorAbs, OnRegion Dom, Format Table, SendToServer "Output/1L2-Error (absolute)", File > tmp5];
+ Print [errorRel, OnRegion Dom, Format Table, SendToServer "Output/2L2-Error (relative)", File > tmp6];
+ }
+ }
+ { Name ProjErrors; NameOfPostProcessing ProjErrors;
+ Operation {
+ tmpA = Sprintf("out/solRef.pos");
+ tmpB = Sprintf("out/solProj.pos");
+ tmpC = Sprintf("out/solDiff.pos");
+ tmp1 = Sprintf("out/errorAbs_%g.dat", FLAG_SIGNAL_BC);
+ tmp2 = Sprintf("out/errorRel_%g.dat", FLAG_SIGNAL_BC);
+ Print [u_ref , OnElementsOf Dom, File tmpA];
+ Print [u_proj, OnElementsOf Dom, File tmpB];
+ Print [u_diff, OnElementsOf Dom, File tmpC];
+ Print [error2[Dom], OnRegion Dom, Format Table, SendToServer "Output/2Error2 (absolute)", StoreInVariable $error2Var];
+ Print [energy2[Dom], OnRegion Dom, Format Table, SendToServer "Output/3Energy2 (absolute)", StoreInVariable $energy2Var];
+ Print [errorAbs, OnRegion Dom, Format Table, SendToServer "Output/1L2-Error (absolute)", File > tmp1];
+ Print [errorRel, OnRegion Dom, Format Table, SendToServer "Output/2L2-Error (relative)", File > tmp2];
+ }
+ }
+}
+
+DefineConstant[
+ R_ = {"NumSol", Name "GetDP/1ResolutionChoices", Visible 1, Choices {"NumNormal", "NumCur", "NumSol", "ProjSol"}},
+ P_ = {"NumSol, Errors", Name "GetDP/2PostOperationChoices", Visible 1, Choices {"NumSol", "Errors", "ProjErrors"}}
+];
diff --git a/HelmholtzHABCwithCorners/padePolyhedron.dat b/HelmholtzHABCwithCorners/padePolyhedron.dat
new file mode 100644
index 0000000000000000000000000000000000000000..5bf4b75ac74937379754a5e6697ca0105fa986e5
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padePolyhedron.dat
@@ -0,0 +1,227 @@
+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;
diff --git a/HelmholtzHABCwithCorners/padePolyhedron.geo b/HelmholtzHABCwithCorners/padePolyhedron.geo
new file mode 100644
index 0000000000000000000000000000000000000000..469ce141dbf2d0105136f86b1486241eb6d56226
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padePolyhedron.geo
@@ -0,0 +1,284 @@
+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;
diff --git a/HelmholtzHABCwithCorners/padePolyhedron.pro b/HelmholtzHABCwithCorners/padePolyhedron.pro
new file mode 100644
index 0000000000000000000000000000000000000000..048368acf3c91e873f3f8bf4f883179384177f03
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padePolyhedron.pro
@@ -0,0 +1,564 @@
+Include "padePolyhedron.dat";
+
+Function {
+If((FLAG_DIM == 3) && (FLAG_SIGNAL == SIGNAL_Scatt))
+If((FAC_Num == 4) || (FAC_Num == 12) || (FAC_Num == 20) || (FAC_Num == 1))
+ incDir = {1., 0., 0.};
+Else
+ incDir = {1./Sqrt[2.], 1./Sqrt[2.], 0.};
+EndIf
+ incPha[] = k[] * (incDir(0)*X[] + incDir(1)*Y[] + incDir(2)*Z[]);
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+ f_inc2[] = Complex[Cos[incPha[]], Sin[incPha[]]];
+ f_ref2[] = AcousticFieldSoftSphere[XYZ[]]{WAVENUMBER, R_SCA, incDir(0), incDir(1), incDir(2)};
+EndIf
+If(FLAG_SIGNAL_BC == SIGNAL_Neumann)
+ df_inc2[] = Complex[0,1] * incPha[]/R[] * Complex[Cos[incPha[]], Sin[incPha[]]];
+ f_ref2[] = AcousticFieldHardSphere[XYZ[]]{WAVENUMBER, R_SCA, incDir(0), incDir(1), incDir(2)};
+EndIf
+EndIf
+}
+
+//==================================================================================================
+// OPTIONS and PARAMETERS
+//==================================================================================================
+
+BND_Neumann = 0;
+BND_Sommerfeld = 1;
+BND_PadeCont = 3;
+BND_PadeDisc = 4;
+CRN_Nothing = 0;
+CRN_Curvature = 2;
+CRN_Curvature2 = 3;
+CRN_SommerfeldOnEdge = 4;
+CRN_ApproxCompOnEdge = 5;
+CRN_SommerfeldAtCorner = 6;
+CRN_ApproxCompAtCorner = 7;
+
+DefineConstant[
+ BND_TYPE = {BND_PadeDisc,
+ Name "Input/5Model/03Boundary condition (faces)", Highlight "Blue",
+ Choices {BND_Neumann = "Homogeneous Neumann",
+ BND_Sommerfeld = "Sommerfeld ABC",
+ BND_PadeDisc = "Pade ABC (disc at corners)",
+ BND_PadeCont = "Pade ABC (cont at corners)"}},
+ CRN_TYPE = {CRN_ApproxCompAtCorner,
+ Name "Input/5Model/04Boundary condition (edges-corners)", Highlight "Red",
+ Visible ((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont)),
+ Choices {CRN_Nothing = "Nothing",
+ CRN_Curvature = "Curvature (kEps)",
+ CRN_Curvature2 = "Curvature (kEps+AddTerms)",
+ CRN_SommerfeldOnEdge = "Sommerfeld on Edges",
+ CRN_ApproxCompOnEdge = "2D-compatibility on Edges",
+ CRN_SommerfeldAtCorner = "Sommerfeld at Corners",
+ CRN_ApproxCompAtCorner = "3D-compatibility at Corners"}}
+];
+
+DefineConstant[
+ nPade = {1, Choices {1, 2, 3, 4},
+ Name "Input/5Model/05Pade: Number of fields",
+ Visible ((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))},
+ thetaPadeInput = {0, Min 0, Step 1, Max 4,
+ Name "Input/5Model/06Pade: Rotation of branch cut",
+ Visible ((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))}
+];
+
+If((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))
+ If(thetaPadeInput == 0)
+ thetaPade = 0;
+ EndIf
+ If(thetaPadeInput == 1)
+ thetaPade = Pi/8;
+ EndIf
+ If(thetaPadeInput == 2)
+ thetaPade = Pi/4;
+ EndIf
+ If(thetaPadeInput == 3)
+ thetaPade = Pi/3;
+ EndIf
+ If(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
+
+Group {
+ Dom = Region[{VOL}];
+ BndSca = Region[{SUR_Scatt}];
+ BndExt = {};
+ For iFac In {1:FAC_Num}
+ BndExt += Region[{SUR~{iFac}}];
+ Face~{iFac} = Region[{SUR~{iFac}}];
+ BndFace~{iFac} = {};
+ EndFor
+ For iEdg In {1:EDG_Num}
+ Edg~{iEdg} = Region[{EDG~{iEdg}}];
+ BndEdg~{iEdg} = {};
+ For iNeigh In {1:EDG_NumFacPerEdg}
+ iFac = EdgNeigh~{iEdg}~{iNeigh};
+ BndFace~{iFac} += Region[{EDG~{iEdg}}];
+ EndFor
+ EndFor
+ For iCrn In {1:CRN_Num}
+ Crn~{iCrn} = Region[{CRN~{iCrn}}];
+ For iNeigh In {1:CRN_NumEdgPerCrn}
+ iEdg = CrnNeigh~{iCrn}~{iNeigh};
+ BndEdg~{iEdg} += Region[{CRN~{iCrn}}];
+ EndFor
+ EndFor
+ DomAll = Region[{Dom,BndSca,BndExt}];
+}
+
+Function {
+ NormalNum[] = VectorField[XYZ[]]{1001};
+ CurvNum[] = ScalarField[XYZ[]]{1002};
+
+If((CRN_TYPE == CRN_Curvature) || (CRN_TYPE == CRN_Curvature2))
+If(FAC_Num == 1)
+ Curv[BndExt] = 1/dimRmean;
+Else
+ Curv[BndExt] = CurvNum[];
+EndIf
+Else
+ Curv[BndExt] = 0.;
+EndIf
+
+ kEps[BndExt] = WAVENUMBER + I[] * 0.4 * WAVENUMBER^(1/3) * (Curv[])^(2/3);
+
+If((BND_TYPE == BND_PadeDisc) || (BND_TYPE == BND_PadeCont))
+ ExpPTheta[] = Complex[Cos[ thetaPade],Sin[ thetaPade]];
+ ExpMTheta[] = Complex[Cos[-thetaPade],Sin[-thetaPade]];
+ ExpPTheta2[] = Complex[Cos[thetaPade/2.],Sin[thetaPade/2.]];
+EndIf
+}
+
+//==================================================================================================
+// FONCTION SPACES with CONSTRAINTS
+//==================================================================================================
+
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+Constraint {
+ { Name DirichletBC; Case {{ Region BndSca; Value -f_inc2[]; }}}
+}
+EndIf
+
+FunctionSpace {
+ { Name H_bndInt; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
+ { Name H_bndExt; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{BndSca}]; Entity NodesOf[All]; }}}
+ { Name H_nx; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
+ { Name H_ny; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
+ { Name H_nz; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
+ { Name H_cur; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
+ { Name H_proj; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{Dom}]; Entity NodesOf[All]; }}}
+ { Name H_num; Type Form0; BasisFunction {{ Name sn; NameOfCoef pn; Function BF_Node; Support Region[{DomAll}]; Entity NodesOf[All]; }}
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+ Constraint {{ NameOfCoef pn; EntityType NodesOf; NameOfConstraint DirichletBC; }}
+EndIf
+ }
+If(BND_TYPE == BND_PadeCont)
+ For m In {1:nPade}
+ { Name H~{m}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
+ EndFor
+EndIf
+If(BND_TYPE == BND_PadeDisc)
+ For iFac In {1:FAC_Num}
+ For m In {1:nPade}
+ { Name H~{iFac}~{m}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{Face~{iFac},BndFace~{iFac}}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+If((CRN_TYPE == CRN_ApproxCompOnEdge) || (CRN_TYPE == CRN_SommerfeldAtCorner) || (CRN_TYPE == CRN_ApproxCompAtCorner))
+ For iEdg In {1:EDG_Num}
+ For m In {1:nPade}
+ For n In {1:nPade}
+ { Name H~{iEdg}~{m}~{n}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{Edg~{iEdg},BndEdg~{iEdg}}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+ EndFor
+If(CRN_TYPE == CRN_ApproxCompAtCorner)
+ For iCrn In {1:CRN_Num}
+ For m In {1:nPade}
+ For n In {1:nPade}
+ For o In {1:nPade}
+ { Name H~{iCrn}~{m}~{n}~{o}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{Crn~{iCrn}}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+ EndFor
+ EndFor
+EndIf
+EndIf
+EndIf
+}
+
+//==================================================================================================
+// FORMULATIONS
+//==================================================================================================
+
+Formulation {
+ { Name NumNormal; Type FemEquation;
+ Quantity {
+ { Name u_nx; Type Local; NameOfSpace H_nx; }
+ { Name u_ny; Type Local; NameOfSpace H_ny; }
+ { Name u_nz; Type Local; NameOfSpace H_nz; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_nx} , {u_nx} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ny} , {u_ny} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_nz} , {u_nz} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompX[Normal[]] , {u_nx} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompY[Normal[]] , {u_ny} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompZ[Normal[]] , {u_nz} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ }
+ }
+
+ { Name NumCur; Type FemEquation;
+ Quantity {
+ { Name u_nx; Type Local; NameOfSpace H_nx; }
+ { Name u_ny; Type Local; NameOfSpace H_ny; }
+ { Name u_nz; Type Local; NameOfSpace H_nz; }
+ { Name u_cur; Type Local; NameOfSpace H_cur; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_nx} , {u_nx} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ny} , {u_ny} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_nz} , {u_nz} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompX[NormalNum[]] , {u_nx} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompY[NormalNum[]] , {u_ny} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompZ[NormalNum[]] , {u_nz} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+
+ Galerkin{ [ Dof{u_cur} , {u_cur} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[0.5,0,0] * Dof{d u_nx} , {u_cur} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[0,0.5,0] * Dof{d u_ny} , {u_cur} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[0,0,0.5] * Dof{d u_nz} , {u_cur} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ }
+ }
+
+ { Name NumSol; Type FemEquation;
+ Quantity {
+ { Name u_num; Type Local; NameOfSpace H_num; }
+If(BND_TYPE == BND_PadeCont)
+ For m In {1:nPade}
+ { Name u~{m}; Type Local; NameOfSpace H~{m}; }
+ EndFor
+EndIf
+If(BND_TYPE == BND_PadeDisc)
+ For iFac In {1:FAC_Num}
+ For m In {1:nPade}
+ { Name u~{iFac}~{m}; Type Local; NameOfSpace H~{iFac}~{m}; }
+ EndFor
+ EndFor
+If((CRN_TYPE == CRN_ApproxCompOnEdge) || (CRN_TYPE == CRN_SommerfeldAtCorner) || (CRN_TYPE == CRN_ApproxCompAtCorner))
+ For iEdg In {1:EDG_Num}
+ For m In {1:nPade}
+ For n In {1:nPade}
+ { Name u~{iEdg}~{m}~{n}; Type Local; NameOfSpace H~{iEdg}~{m}~{n}; }
+ EndFor
+ EndFor
+ EndFor
+If(CRN_TYPE == CRN_ApproxCompAtCorner)
+ For iCrn In {1:CRN_Num}
+ For m In {1:nPade}
+ For n In {1:nPade}
+ For o In {1:nPade}
+ { Name u~{iCrn}~{m}~{n}~{o}; Type Local; NameOfSpace H~{iCrn}~{m}~{n}~{o}; }
+ EndFor
+ EndFor
+ EndFor
+ EndFor
+EndIf
+EndIf
+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_inc2[] , {u_num} ]; In BndSca; Jacobian JSur; 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
+
+// HABC (continuity at corners)
+If(BND_TYPE == BND_PadeCont)
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+ If(CRN_TYPE == CRN_Curvature2)
+ Galerkin { [ Curv[] * Dof{u_num} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+ Galerkin { [ -Curv[]/(2*k[]*k[]) * Dof{d u_num} , {d u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+ EndIf
+ For m In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{m}} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u_num} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u~{m}} , {d u~{m}} ]; In BndExt; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{m}+ExpMTheta[]) * Dof{u~{m}} , {u~{m}} ]; In BndExt; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{m}+1) * Dof{u_num} , {u~{m}} ]; In BndExt; Jacobian JSur; Integration I1; }
+ EndFor
+EndIf // End HABC Cont
+
+// HABC (discontinuity at corners)
+If(BND_TYPE == BND_PadeDisc)
+ For iFac In {1:FAC_Num}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In Face~{iFac}; Jacobian JSur; Integration I1; }
+ For m In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{iFac}~{m}} , {u_num} ]; In Face~{iFac}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u_num} , {u_num} ]; In Face~{iFac}; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u~{iFac}~{m}} , {d u~{iFac}~{m}} ]; In Face~{iFac}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+ExpMTheta[]) * Dof{u~{iFac}~{m}} , {u~{iFac}~{m}} ]; In Face~{iFac}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+1) * Dof{u_num} , {u~{iFac}~{m}} ]; In Face~{iFac}; Jacobian JSur; Integration I1; }
+
+ If(CRN_TYPE == CRN_SommerfeldOnEdge)
+ Galerkin { [ -I[]*k[] * Dof{u~{iFac}~{m}} , {u~{iFac}~{m}} ]; In BndFace~{iFac}; Jacobian JLin; Integration I1; }
+ EndIf
+ EndFor
+ EndFor
+
+If((CRN_TYPE == CRN_ApproxCompOnEdge) || (CRN_TYPE == CRN_SommerfeldAtCorner) || (CRN_TYPE == CRN_ApproxCompAtCorner)) /////
+ For iEdg In {1:EDG_Num}
+ iFace1 = EdgNeigh~{iEdg}~{1};
+ iFace2 = EdgNeigh~{iEdg}~{2};
+ For m In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iFace1}~{m}} , {u~{iFace1}~{m}} ]; In Edg~{iEdg}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iFace2}~{m}} , {u~{iFace2}~{m}} ]; In Edg~{iEdg}; Jacobian JLin; Integration I1; }
+
+ For n In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{iFace1}~{m}} , {u~{iFace1}~{m}} ]; In Edg~{iEdg}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{iFace2}~{m}} , {u~{iFace2}~{m}} ]; In Edg~{iEdg}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{iEdg}~{m}~{n}} , {u~{iFace1}~{m}} ]; In Edg~{iEdg}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{iEdg}~{n}~{m}} , {u~{iFace2}~{m}} ]; In Edg~{iEdg}; Jacobian JLin; Integration I1; }
+
+ If((CRN_TYPE == CRN_SommerfeldAtCorner) || (CRN_TYPE == CRN_ApproxCompAtCorner))
+ Galerkin { [ Dof{d u~{iEdg}~{m}~{n}} , {d u~{iEdg}~{m}~{n}} ]; In Edg~{iEdg}; Jacobian JLin; Integration I1; }
+ EndIf
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+cPade~{n}+ExpMTheta[]) * Dof{u~{iEdg}~{m}~{n}} , {u~{iEdg}~{m}~{n}} ]; In Edg~{iEdg}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{n}+1) * Dof{u~{iFace1}~{m}} , {u~{iEdg}~{m}~{n}} ]; In Edg~{iEdg}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -k[]^2*ExpPTheta[]*(cPade~{m}+1) * Dof{u~{iFace2}~{n}} , {u~{iEdg}~{m}~{n}} ]; In Edg~{iEdg}; Jacobian JLin; Integration I1; }
+
+ If(CRN_TYPE == CRN_SommerfeldAtCorner)
+ Galerkin { [ -I[]*k[] * Dof{u~{iEdg}~{m}~{n}} , {u~{iEdg}~{m}~{n}} ]; In BndEdg~{iEdg}; Jacobian JVol; Integration I1; }
+ EndIf
+ EndFor
+ EndFor
+ EndFor
+EndIf
+
+If(CRN_TYPE == CRN_ApproxCompAtCorner)
+ For iCrn In {1:CRN_Num}
+ iEdg1 = CrnNeigh~{iCrn}~{1}; // 1-2 m-n
+ iEdg2 = CrnNeigh~{iCrn}~{2}; // 1-3 m-o
+ iEdg3 = CrnNeigh~{iCrn}~{3}; // 2-3 n-o
+ For m In {1:nPade}
+ For n In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iEdg1}~{m}~{n}} , {u~{iEdg1}~{m}~{n}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iEdg2}~{m}~{n}} , {u~{iEdg2}~{m}~{n}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iEdg3}~{m}~{n}} , {u~{iEdg3}~{m}~{n}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+
+ For o In {1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{o} * Dof{u~{iCrn}~{m}~{n}~{o}} , {u~{iEdg1}~{m}~{n}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{o} * Dof{u~{iEdg1}~{m}~{n}} , {u~{iEdg1}~{m}~{n}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{iCrn}~{m}~{n}~{o}} , {u~{iEdg2}~{m}~{o}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{n} * Dof{u~{iEdg2}~{m}~{o}} , {u~{iEdg2}~{m}~{o}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{iCrn}~{m}~{n}~{o}} , {u~{iEdg3}~{n}~{o}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{m} * Dof{u~{iEdg3}~{n}~{o}} , {u~{iEdg3}~{n}~{o}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+
+ Galerkin { [ (cPade~{m}+cPade~{n}+cPade~{o}+ExpMTheta[]) * Dof{u~{iCrn}~{m}~{n}~{o}} , {u~{iCrn}~{m}~{n}~{o}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+ Galerkin { [ (cPade~{o}+1) * Dof{u~{iEdg1}~{m}~{n}} , {u~{iCrn}~{m}~{n}~{o}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+ Galerkin { [ (cPade~{n}+1) * Dof{u~{iEdg2}~{m}~{o}} , {u~{iCrn}~{m}~{n}~{o}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+ Galerkin { [ (cPade~{m}+1) * Dof{u~{iEdg3}~{n}~{o}} , {u~{iCrn}~{m}~{n}~{o}} ]; In Crn~{iCrn}; Jacobian JVol; Integration I1; }
+ EndFor
+ EndFor
+ EndFor
+ EndFor
+EndIf
+
+EndIf // End HABC Cont
+
+ }
+ }
+
+ { Name ProjSol; Type FemEquation;
+ Quantity {
+ { Name u_proj; Type Local; NameOfSpace H_proj; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_proj} , {u_proj} ]; In Dom; Jacobian JVol; Integration I1; }
+ Galerkin{ [ -f_ref2[] , {u_proj} ]; In Dom; Jacobian JVol; Integration I1; }
+ }
+ }
+
+ { Name NumBnd; Type FemEquation;
+ Quantity {
+ { Name u_int; Type Local; NameOfSpace H_bndInt; }
+ { Name u_ext; Type Local; NameOfSpace H_bndExt; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_int} , {u_int} ]; In Region[{BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ext} , {u_ext} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -1 , {u_int} ]; In Region[{BndSca}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -1 , {u_ext} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ }
+ }
+}
+
+//==================================================================================================
+// RESOLUTION
+//==================================================================================================
+
+Resolution {
+ { Name NumBnd;
+ System {{ Name A; NameOfFormulation NumBnd; Type Real; }}
+ Operation { Generate[A]; Solve[A]; SaveSolution[A]; }
+ }
+ { Name NumNormal;
+ System {
+ { Name A; NameOfFormulation NumNormal; Type Real; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A];
+ }
+ }
+ { Name NumCur;
+ System {
+ { Name A; NameOfFormulation NumNormal; Type Real; }
+ { Name B; NameOfFormulation NumCur; Type Real; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A]; PostOperation[NumNormal];
+ Generate[B]; Solve[B]; SaveSolution[B];
+ }
+ }
+ { Name NumSol;
+ System {
+If(((CRN_TYPE == CRN_Curvature) || (CRN_TYPE == CRN_Curvature2)) && (FAC_Num > 1))
+ { Name A; NameOfFormulation NumNormal; Type Real; NameOfMesh "mainCurv.msh"; }
+ { Name B; NameOfFormulation NumCur; Type Real; NameOfMesh "mainCurv.msh"; }
+EndIf
+ { Name C; NameOfFormulation NumSol; Type Complex; }
+ }
+ Operation {
+If(((CRN_TYPE == CRN_Curvature) || (CRN_TYPE == CRN_Curvature2)) && (FAC_Num > 1))
+ Generate[A]; Solve[A]; SaveSolution[A]; PostOperation[NumNormal];
+ Generate[B]; Solve[B]; SaveSolution[B]; PostOperation[NumCur];
+EndIf
+ Generate[C]; Solve[C]; SaveSolution[C];
+ }
+ }
+ { Name ProjSol;
+ System {{ Name A; NameOfFormulation ProjSol; Type Complex; }}
+ Operation { Generate[A]; Solve[A]; SaveSolution[A]; }
+ }
+}
+
+//==================================================================================================
+// POSTPRO / POSTOP
+//==================================================================================================
+
+PostProcessing {
+ { Name NumNormal; NameOfFormulation NumNormal;
+ Quantity {
+ { Name u_normal; Value { Local { [ Vector[{u_nx},{u_ny},{u_nz}] / Sqrt[{u_nx}*{u_nx}+{u_ny}*{u_ny}+{u_nz}*{u_nz}] ]; In Region[{BndExt}]; Jacobian JSur; }}}
+ }
+ }
+ { Name NumCur; NameOfFormulation NumCur;
+ Quantity {
+ { Name u_cur; Value { Local { [ {u_cur} ]; In Region[{BndExt}]; Jacobian JSur; }}}
+ }
+ }
+ { Name NumSol; NameOfFormulation NumSol;
+ Quantity {
+ { Name u_ref~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ f_ref2[] ]; In Dom; Jacobian JVol; }}}
+ { Name u_num~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ {u_num} ]; In Dom; Jacobian JVol; }}}
+ { Name u_err~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ f_ref2[]-{u_num}]; In Dom; Jacobian JVol; }}}
+ }
+ }
+ { Name Errors; NameOfFormulation NumSol;
+ Quantity {
+ { Name error2; Value { Integral { [ Norm[f_ref2[]-{u_num}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[f_ref2[]]^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/$energy2Var]] ; In Dom; Jacobian JVol; }}}
+ }
+ }
+ { Name ProjError; NameOfFormulation ProjSol;
+ Quantity {
+ { Name error2; Value { Integral { [ Norm[f_ref2[]-{u_proj}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
+ { Name energy2; Value { Integral { [ Norm[f_ref2[]]^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/$energy2Var]] ; In Dom; Jacobian JVol; }}}
+ }
+ }
+ { Name NumBnd; NameOfFormulation NumBnd;
+ Quantity {
+ { Name u_int; Value { Local { [ {u_int} ]; In BndSca; Jacobian JSur; }}}
+ { Name u_ext; Value { Local { [ {u_ext} ]; In BndExt; Jacobian JSur; }}}
+ }
+ }
+}
+
+PostOperation{
+ { Name NumNormal; NameOfPostProcessing NumNormal;
+ Operation {
+ Print [u_normal, OnElementsOf Region[{BndExt}], StoreInField (1001), File "out/u_normal.pos"];
+ }
+ }
+ { Name NumCur; NameOfPostProcessing NumCur;
+ Operation {
+ Print [u_cur, OnElementsOf Region[{BndExt}], StoreInField (1002), File "out/u_cur.pos"];
+ }
+ }
+ { Name NumSol; NameOfPostProcessing NumSol;
+ Operation {
+ tmp1 = StrCat(DIR, Sprintf("solRef_poly_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, FAC_Num));
+ tmp2 = StrCat(DIR, Sprintf("solNum_poly_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, FAC_Num));
+ tmp3 = StrCat(DIR, Sprintf("solErr_poly_%g_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, FAC_Num));
+ Print [u_ref~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Dom, File tmp1];
+ Print [u_num~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Dom, File tmp2];
+ Print [u_err~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Dom, File tmp3];
+ }
+ }
+ { Name Errors; NameOfPostProcessing Errors;
+ Operation {
+ tmp4 = StrCat(DIR, Sprintf("errorAbs_poly_%g_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, FAC_Num));
+ tmp5 = StrCat(DIR, Sprintf("errorRel_poly_%g_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput, FAC_Num));
+ Print [error2[Dom], OnRegion Dom, Format Table, StoreInVariable $error2Var];
+ Print [energy2[Dom], OnRegion Dom, Format Table, StoreInVariable $energy2Var];
+ Print [errorAbs, OnRegion Dom, Format Table, SendToServer "Output/1L2-Error (absolute)", File > tmp4];
+ Print [errorRel, OnRegion Dom, Format Table, SendToServer "Output/2L2-Error (relative)", File > tmp5];
+ }
+ }
+ { Name ProjError; NameOfPostProcessing ProjError;
+ Operation {
+ tmp6 = StrCat(DIR, Sprintf("errorAbs_%g_%g_%g.dat", FLAG_SIGNAL_BC, FAC_Num, N_LAMBDA));
+ tmp7 = StrCat(DIR, Sprintf("errorRel_%g_%g_%g.dat", FLAG_SIGNAL_BC, FAC_Num, N_LAMBDA));
+ Print [error2[Dom], OnRegion Dom, Format Table, StoreInVariable $error2Var];
+ Print [energy2[Dom], OnRegion Dom, Format Table, StoreInVariable $energy2Var];
+ Print [errorAbs, OnRegion Dom, Format Table, SendToServer "Output/1L2-Error (absolute)", File > tmp6];
+ Print [errorRel, OnRegion Dom, Format Table, SendToServer "Output/2L2-Error (relative)", File > tmp7];
+ }
+ }
+ { Name NumBnd; NameOfPostProcessing NumBnd;
+ Operation {
+ tmp1 = StrCat(DIR, Sprintf("bndInt_poly_%g_%g.pos", FLAG_SIGNAL_BC, FAC_Num));
+ tmp2 = StrCat(DIR, Sprintf("bndExt_poly_%g_%g.pos", FLAG_SIGNAL_BC, FAC_Num));
+ Print [u_int, OnElementsOf BndSca, File tmp1];
+ Print [u_ext, OnElementsOf BndExt, File tmp2];
+ }
+ }
+}
+
+DefineConstant[
+ R_ = {"NumSol", Name "GetDP/1ResolutionChoices", Visible 1, Choices {"NumNormal", "NumCur", "NumSol", "ProjSol", "NumBnd"}},
+ P_ = {"NumSol, Errors", Name "GetDP/2PostOperationChoices", Visible 1, Choices {"NumNormal", "NumCur", "NumSol", "Errors", "ProjError", "NumBnd"}}
+];
diff --git a/HelmholtzHABCwithCorners/padeSquare.dat b/HelmholtzHABCwithCorners/padeSquare.dat
new file mode 100644
index 0000000000000000000000000000000000000000..da871f78bce344a4332a5910a91427c9c16bdac0
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padeSquare.dat
@@ -0,0 +1,53 @@
+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;
diff --git a/HelmholtzHABCwithCorners/padeSquare.geo b/HelmholtzHABCwithCorners/padeSquare.geo
new file mode 100644
index 0000000000000000000000000000000000000000..08ab0bfdce12ae3af88a5bf6147be53ffd351a6a
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padeSquare.geo
@@ -0,0 +1,71 @@
+Include "padeSquare.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(4) = {X~{4}, Y~{4}, 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};
+
+Line(1) = {4, 1};
+Line(2) = {1, 2};
+Line(3) = {2, 3};
+Line(4) = {3, 4};
+
+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)
+ out0[] = Extrude {-Lpml,0,0} { Line{1}; Layers{Npml}; Recombine; };
+ Printf("==== %g ====\n", out0[0]);
+ Printf("==== %g ====\n", out0[1]);
+ Printf("==== %g ====\n", out0[2]);
+ Printf("==== %g ====\n", out0[3]);
+ out1[] = Extrude { Lpml,0,0} { Line{3}; Layers{Npml}; Recombine; };
+ Printf("==== %g ====\n", out1[0]);
+ Printf("==== %g ====\n", out1[1]);
+ Printf("==== %g ====\n", out1[2]);
+ Printf("==== %g ====\n", out1[3]);
+ out2[] = Extrude {0,-Lpml,0} { Line{11,2,14}; Layers{Npml}; Recombine; };
+ Printf("==== %g ====\n", out2[0]);
+ Printf("==== %g ====\n", out2[1]);
+ Printf("==== %g ====\n", out2[2]);
+ Printf("==== %g ====\n", out2[3]);
+ out3[] = Extrude {0, Lpml,0} { Line{10,4,15}; Layers{Npml}; Recombine; };
+ Printf("==== %g ====\n", out3[0]);
+ Printf("==== %g ====\n", out3[1]);
+ Printf("==== %g ====\n", out3[2]);
+ Printf("==== %g ====\n", out3[3]);
+EndIf
+
+Physical Point(CRN_1) = {1};
+Physical Point(CRN_2) = {2};
+Physical Point(CRN_3) = {3};
+Physical Point(CRN_4) = {4};
+
+Physical Line(BND_1) = {1}; // Left
+Physical Line(BND_2) = {2}; // Down
+Physical Line(BND_3) = {3}; // Right
+Physical Line(BND_4) = {4}; // Top
+Physical Line(BND_Scatt) = {5,6,7,8};
+
+Physical Surface(DOM) = {1};
+
+If(BND_TYPE == BND_PML)
+ Physical Surface(DOM_PML) = {12,16,20,24,28,32,36,40};
+If(Npml == 0)
+ Physical Line(BND_PmlExt) = {5,6,7,8};
+Else
+ Physical Line(BND_PmlExt) = {17,21,25,27,13,39,37,33,29,31,9,19};
+EndIf
+EndIf
diff --git a/HelmholtzHABCwithCorners/padeSquare.pro b/HelmholtzHABCwithCorners/padeSquare.pro
new file mode 100644
index 0000000000000000000000000000000000000000..3e900e5976f5163fbad4182b6182fa337fe92c48
--- /dev/null
+++ b/HelmholtzHABCwithCorners/padeSquare.pro
@@ -0,0 +1,565 @@
+Include "padeSquare.dat";
+
+//==================================================================================================
+// OPTIONS and PARAMETERS
+//==================================================================================================
+
+CRN_Regularization = 0;
+CRN_Compatibility = 1;
+KEPS_Nothing = 0;
+KEPS_Analytic = 1;
+KEPS_Numeric = 2;
+
+DefineConstant[
+ CRN_TYPE = {CRN_Compatibility,
+ Name "Input/5Model/03Boundary condition (corners)", Highlight "Red",
+ Visible ((BND_TYPE == BND_Second) || (BND_TYPE == BND_Pade) || (BND_TYPE == BND_CRBC)),
+ Choices {CRN_Regularization = "Regularization",
+ CRN_Compatibility = "Compatibility"}},
+ nPade = {4, Choices {0, 1, 2, 3, 4, 5, 6},
+ Name "Input/5Model/05Pade: Number of fields",
+ Visible ((BND_TYPE == BND_Pade) || (BND_TYPE == BND_CRBC))},
+ thetaPadeInput = {3, Choices {0, 1, 2, 3, 4},
+ Name "Input/5Model/06Pade: Rotation of branch cut",
+ Visible ((BND_TYPE == BND_Pade) || (BND_TYPE == BND_CRBC))},
+ KEPS_TYPE = {KEPS_Numeric,
+ Name "Input/5Model/07Curvature for regularization", Highlight "Red",
+ Visible (((BND_TYPE == BND_Pade) || (BND_TYPE == BND_CRBC)) && (CRN_TYPE == CRN_Regularization)),
+ Choices {KEPS_Nothing = "Nothing",
+ KEPS_Analytic = "Analytic formula",
+ KEPS_Numeric = "Numerical curvature"}}
+];
+
+If((BND_TYPE == BND_Pade) || (BND_TYPE == BND_CRBC))
+ If(thetaPadeInput == 0)
+ thetaPade = 0;
+ EndIf
+ If(thetaPadeInput == 1)
+ thetaPade = Pi/8;
+ EndIf
+ If(thetaPadeInput == 2)
+ thetaPade = Pi/4;
+ EndIf
+ If(thetaPadeInput == 3)
+ thetaPade = Pi/3;
+ EndIf
+ If(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}];
+If(BND_TYPE == BND_PML)
+ DomPml = Region[{DOM_PML}];
+ BndPml = Region[{BND_PmlExt}];
+ DomPmlAll = Region[{DomPml,BndPml}];
+Else
+ DomPml = Region[{}];
+ BndPml = Region[{}];
+ DomPmlAll = Region[{}];
+EndIf
+ BndSca = Region[{BND_Scatt}];
+ For iEdg In{1:4}
+ iCrn1 = (iEdg == 1) ? 4 : iEdg-1;
+ iCrn2 = iEdg;
+ Edg~{iEdg} = Region[{BND~{iEdg}}];
+ EdgClo~{iEdg} = Region[{BND~{iEdg},CRN~{iCrn1},CRN~{iCrn2}}];
+ EndFor
+ For iCrn In{1:4}
+ Crn~{iCrn} = Region[{CRN~{iCrn}}];
+ EndFor
+
+ CrnAll = Region[{CRN~{1},CRN~{2},CRN~{3},CRN~{4}}];
+ EdgAll = Region[{BND~{1},BND~{2},BND~{3},BND~{4}}];
+ DomAll = Region[{Dom,BndSca,EdgAll,CrnAll}];
+ BndExt = Region[{EdgAll}];
+}
+
+Function {
+
+ For iEdg In{1:4}
+ NormalGeo[Edg~{iEdg}] = Vector[Cos[(iEdg+1)*Pi/2], Sin[(iEdg+1)*Pi/2], 0.];
+ EndFor
+
+ NormalNum[] = VectorField[XYZ[]]{1001};
+ CurvNum[] = ScalarField[XYZ[]]{1002};
+
+If(((BND_TYPE == BND_Pade) || (BND_TYPE == BND_CRBC)) && (CRN_TYPE == CRN_Regularization) && (KEPS_TYPE == KEPS_Analytic))
+ KAPPA = (2)^(1/2)/LC;
+ dist~{1}[] = Sqrt[(X[]-X~{1})^2 + (Y[]-Y~{1})^2];
+ dist~{2}[] = Sqrt[(X[]-X~{2})^2 + (Y[]-Y~{2})^2];
+ dist~{3}[] = Sqrt[(X[]-X~{3})^2 + (Y[]-Y~{3})^2];
+ dist~{4}[] = Sqrt[(X[]-X~{4})^2 + (Y[]-Y~{4})^2];
+ Curv[EdgAll] = ((dist~{1}[] < LC)||(dist~{2}[] < LC)||(dist~{3}[] < LC)||(dist~{4}[] < LC)) ? KAPPA : 0;
+ kEps[] = WAVENUMBER + I[] * 0.39 * WAVENUMBER^(1/3) * (Curv[])^(2/3);
+ElseIf(((BND_TYPE == BND_Pade) || (BND_TYPE == BND_CRBC)) && (CRN_TYPE == CRN_Regularization) && (KEPS_TYPE == KEPS_Numeric))
+ kEps[] = WAVENUMBER + I[] * 0.39 * WAVENUMBER^(1/3) * (CurvNum[])^(2/3);
+Else
+ Curv[EdgAll] = 0.;
+ kEps[EdgAll] = WAVENUMBER;
+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_CRBC)
+ For n In{0:nPade}
+ Alpha~{n}[] = Complex[Cos[thetaPade/2.], Sin[thetaPade/2.]];
+ EndFor
+EndIf
+
+If(BND_TYPE == BND_PML)
+ xLoc[DomPml] = Fabs[X[]+X_SCA-dimL/2]-dimL/2;
+ yLoc[DomPml] = Fabs[Y[]+Y_SCA-dimL/2]-dimL/2;
+ //absFuncX[DomPml] = (xLoc[]>=0) ? sigmaPml*(Lpml-xLoc[])*(Lpml-xLoc[]) : 0;
+ //absFuncY[DomPml] = (yLoc[]>=0) ? sigmaPml*(Lpml-yLoc[])*(Lpml-yLoc[]) : 0;
+ absFuncX[DomPml] = (xLoc[]>=0) ? 1/(Lpml-xLoc[]) : 0;
+ absFuncY[DomPml] = (yLoc[]>=0) ? 1/(Lpml-yLoc[]) : 0;
+ 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
+ hx[DomPml] = 1 + rot[] * absFuncX[]/k[];
+ hy[DomPml] = 1 + rot[] * absFuncY[]/k[];
+ pmlScal[DomPml] = hx[]*hy[];
+ pmlTens[DomPml] = TensorDiag[hy[]/hx[], hx[]/hy[], 0];
+EndIf
+}
+
+//==================================================================================================
+// FONCTION SPACES with CONSTRAINTS
+//==================================================================================================
+
+If((FLAG_SIGNAL_BC == SIGNAL_Dirichlet) || (BND_TYPE == BND_PML))
+Constraint {
+ { Name DirichletBC; Case {
+If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
+ { Region BndSca; Value -f_inc[]; }
+EndIf
+ }}
+}
+EndIf
+
+FunctionSpace {
+ { Name H_nx; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
+ { Name H_ny; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
+ { Name H_cur; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
+ { 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,DomPmlAll}]; 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) || (BND_TYPE == BND_CRBC))
+If(CRN_TYPE == CRN_Regularization)
+ For i In {1:nPade}
+ { Name H~{i}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{EdgAll}]; Entity NodesOf[All]; }}}
+ EndFor
+EndIf
+If(CRN_TYPE == CRN_Compatibility)
+ For iEdg In{1:4}
+ For i In {1:nPade}
+ { Name H~{iEdg}~{i}; Type Form0;BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{EdgClo~{iEdg}}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+ For iCrn In{1:4}
+ For i In {1:nPade}
+ For j In {1:nPade}
+ { Name H~{iCrn}~{i}~{j}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{Crn~{iCrn}}]; Entity NodesOf[All]; }}}
+ EndFor
+ EndFor
+ EndFor
+EndIf
+EndIf
+}
+
+//==================================================================================================
+// FORMULATIONS
+//==================================================================================================
+
+Formulation {
+ { Name NumNormal; Type FemEquation;
+ Quantity {
+ { Name u_nx; Type Local; NameOfSpace H_nx; }
+ { Name u_ny; Type Local; NameOfSpace H_ny; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_nx} , {u_nx} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ny} , {u_ny} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompX[NormalGeo[]] , {u_nx} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompY[NormalGeo[]] , {u_ny} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ }
+ }
+
+ { Name NumCur; Type FemEquation;
+ Quantity {
+ { Name u_nx; Type Local; NameOfSpace H_nx; }
+ { Name u_ny; Type Local; NameOfSpace H_ny; }
+ { Name u_cur; Type Local; NameOfSpace H_cur; }
+ }
+ Equation {
+ Galerkin{ [ Dof{u_nx} , {u_nx} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ Dof{u_ny} , {u_ny} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompX[NormalNum[]] , {u_nx} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -CompY[NormalNum[]] , {u_ny} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+
+ Galerkin{ [ Dof{u_cur} , {u_cur} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[1,0,0] * Dof{d u_nx} , {u_cur} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ Galerkin{ [ -Vector[0,1,0] * Dof{d u_ny} , {u_cur} ]; In Region[{BndExt}]; Jacobian JSur; Integration I1; }
+ }
+ }
+
+ { Name NumSol; Type FemEquation;
+ Quantity {
+ { Name u_num; Type Local; NameOfSpace H_num; }
+If((BND_TYPE == BND_Pade) || (BND_TYPE == BND_CRBC))
+If(CRN_TYPE == CRN_Regularization)
+ For i In {1:nPade}
+ { Name u~{i}; Type Local; NameOfSpace H~{i}; }
+ EndFor
+EndIf
+If(CRN_TYPE == CRN_Compatibility)
+ For iEdg In{1:4}
+ For i In {1:nPade}
+ { Name u~{iEdg}~{i}; Type Local; NameOfSpace H~{iEdg}~{i}; }
+ EndFor
+ EndFor
+ For iCrn In{1:4}
+ For i In {1:nPade}
+ For j In {1:nPade}
+ { Name u~{iCrn}~{i}~{j}; Type Local; NameOfSpace H~{iCrn}~{i}~{j}; }
+ EndFor
+ EndFor
+ EndFor
+EndIf
+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; }
+// Galerkin{ [ -I[]*k[]*Dof{u_num} , {u_num} ]; In BndPml; Jacobian JSur; Integration I1; }
+EndIf
+
+// Sommerfeld ABC
+
+If(BND_TYPE == BND_Sommerfeld)
+ Galerkin{ [ -I[]*k[]*Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+EndIf
+
+// Second-order ABC
+
+If(BND_TYPE == BND_Second)
+ Galerkin { [ - I[]*k[] * Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ - 1/(2*I[]*kEps[]) * Dof{d u_num} , {d u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+If(CRN_TYPE == CRN_Compatibility)
+ Galerkin { [ 3./4. * Dof{u_num} , {u_num} ]; In CrnAll; Jacobian JLin; Integration I1; }
+EndIf
+EndIf
+
+// HABC (continuity at corners)
+
+If((BND_TYPE == BND_Pade) && (CRN_TYPE == CRN_Regularization))
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u~{i}} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u_num} , {u_num} ]; In EdgAll; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u~{i}} , {d u~{i}} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*(ExpPTheta[]*cPade~{i}+1) * Dof{u~{i}} , {u~{i}} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{i}+1) * Dof{u_num} , {u~{i}} ]; In EdgAll; Jacobian JSur; Integration I1; }
+ EndFor
+EndIf
+
+// HABC (discontinuity at corners)
+
+If((BND_TYPE == BND_Pade) && (CRN_TYPE == CRN_Compatibility))
+
+ For iEdg In{1:4}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u~{iEdg}~{i}} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u_num} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+
+ Galerkin { [ Dof{d u~{iEdg}~{i}} , {d u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*(ExpPTheta[]*cPade~{i}+1) * Dof{u~{iEdg}~{i}} , {u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{i}+1) * Dof{u_num} , {u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ EndFor
+ EndFor
+
+ For iCrn In{1:4}
+ iEdg1 = iCrn;
+ iEdg2 = (iCrn == 4) ? 1 : iCrn+1;
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ For j In{1:nPade}
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iCrn}~{i}~{j}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{j} * Dof{u~{iCrn}~{j}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+
+ Galerkin { [ (cPade~{i}+cPade~{j}+ExpMTheta[]) * Dof{u~{iCrn}~{i}~{j}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ (cPade~{j}+1) * Dof{u~{iEdg1}~{i}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ (cPade~{i}+1) * Dof{u~{iEdg2}~{j}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ EndFor
+
+// For j In{1:nPade} // Alternative (equivalent) form
+// Galerkin { [ -I[]*k[]*ExpPTheta2[] * coefA~{i}~{j}[] * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+// Galerkin { [ -I[]*k[]*ExpPTheta2[] * coefA~{i}~{j}[] * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+// Galerkin { [ -I[]*k[]*ExpPTheta2[] * coefB~{i}~{j}[] * Dof{u~{iEdg2}~{j}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+// Galerkin { [ -I[]*k[]*ExpPTheta2[] * coefB~{i}~{j}[] * Dof{u~{iEdg1}~{j}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+// EndFor
+
+ EndFor
+ EndFor
+
+EndIf
+
+// CRBC
+
+If((BND_TYPE == BND_CRBC) && (CRN_TYPE == CRN_Compatibility))
+
+ For iEdg In{1:4}
+ Galerkin { [ -I[]*k[] * Alpha~{nPade}[] * Dof{u_num} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ For i1 In{1:nPade}
+ Galerkin { [ -I[]*k[] * (Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u_num} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ For i2 In{i1:nPade}
+ Galerkin { [ -I[]*k[] * (Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u~{iEdg}~{i2}} , {u_num} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ EndFor
+ EndFor
+
+ For i In{1:nPade}
+ Galerkin { [ Dof{d u~{iEdg}~{i}} , {d u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ Galerkin { [ -kEps[]^2 * (1-Alpha~{i}[]^2) * Dof{u~{iEdg}~{i}} , {u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ For i1 In{1:i}
+ Galerkin { [ -2*kEps[]^2 * Alpha~{i}[]*(Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u_num} , {u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ For i2 In{i1:nPade}
+ Galerkin { [ -2*kEps[]^2 * Alpha~{i}[]*(Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u~{iEdg}~{i2}} , {u~{iEdg}~{i}} ]; In Edg~{iEdg}; Jacobian JSur; Integration I1; }
+ EndFor
+ EndFor
+ EndFor
+ EndFor
+
+ For iCrn In{1:4}
+ iEdg1 = iCrn;
+ iEdg2 = (iCrn == 4) ? 1 : iCrn+1;
+
+ For i In{1:nPade}
+ Galerkin { [ -I[]*k[] * Alpha~{nPade}[] * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[] * Alpha~{nPade}[] * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ For i1 In{1:nPade}
+ Galerkin { [ -I[]*k[] * (Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u~{iEdg1}~{i}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[] * (Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u~{iEdg2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ For i2 In{i1:nPade}
+ Galerkin { [ -I[]*k[] * (Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u~{iCrn}~{i}~{i2}} , {u~{iEdg1}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -I[]*k[] * (Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u~{iCrn}~{i2}~{i}} , {u~{iEdg2}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ EndFor
+ EndFor
+ EndFor
+
+ For i In{1:nPade}
+ For j In{1:nPade}
+ Galerkin { [ Dof{u~{iCrn}~{i}~{j}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -(1-Alpha~{i}[]^2) * Dof{u~{iCrn}~{i}~{j}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -(1-Alpha~{i}[]^2) * Dof{u~{iCrn}~{j}~{i}} , {u~{iCrn}~{j}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ For i1 In{1:j}
+ Galerkin { [ -2 * Alpha~{i}[]*(Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u~{iEdg1}~{i}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -2 * Alpha~{i}[]*(Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u~{iEdg2}~{i}} , {u~{iCrn}~{j}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ For i2 In{i1:nPade}
+ Galerkin { [ -2 * Alpha~{i}[]*(Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u~{iCrn}~{i}~{i2}} , {u~{iCrn}~{i}~{j}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ Galerkin { [ -2 * Alpha~{i}[]*(Alpha~{i1}[]+Alpha~{i1-1}[]) * Dof{u~{iCrn}~{i2}~{i}} , {u~{iCrn}~{j}~{i}} ]; In Crn~{iCrn}; Jacobian JLin; Integration I1; }
+ EndFor
+ EndFor
+
+ EndFor
+ EndFor
+
+ 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 NumNormal;
+ System {
+ { Name A; NameOfFormulation NumNormal; Type Real; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A];
+ }
+ }
+ { Name NumCur;
+ System {
+ { Name A; NameOfFormulation NumNormal; Type Real; }
+ { Name B; NameOfFormulation NumCur; Type Real; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A]; PostOperation[NumNormal];
+ Generate[B]; Solve[B]; SaveSolution[B];
+ }
+ }
+ { Name NumSol;
+ System {
+If(((BND_TYPE == BND_Pade) || (BND_TYPE == BND_CRBC)) && (CRN_TYPE == CRN_Regularization) && (KEPS_TYPE == KEPS_Numeric))
+ { Name A; NameOfFormulation NumNormal; Type Real; }
+ { Name B; NameOfFormulation NumCur; Type Real; }
+EndIf
+ { Name C; NameOfFormulation NumSol; Type Complex; }
+ }
+ Operation {
+If(((BND_TYPE == BND_Pade) || (BND_TYPE == BND_CRBC)) && (CRN_TYPE == CRN_Regularization) && (KEPS_TYPE == KEPS_Numeric))
+ Generate[A]; Solve[A]; SaveSolution[A]; PostOperation[NumNormal];
+ Generate[B]; Solve[B]; SaveSolution[B]; PostOperation[NumCur];
+EndIf
+ Generate[C]; Solve[C]; SaveSolution[C];
+ }
+ }
+ { Name ProjSol;
+ System {
+ { Name A; NameOfFormulation ProjSol; Type Complex; }
+ }
+ Operation {
+ Generate[A]; Solve[A]; SaveSolution[A];
+ }
+ }
+}
+
+//==================================================================================================
+// POSTPRO / POSTOP
+//==================================================================================================
+
+PostProcessing {
+ { Name NumNormal; NameOfFormulation NumNormal;
+ Quantity {
+ { Name u_normal; Value { Local { [ Vector[{u_nx},{u_ny},0] / Sqrt[{u_nx}*{u_nx}+{u_ny}*{u_ny}] ]; In Region[{BndExt}]; Jacobian JSur; }}}
+ }
+ }
+ { Name NumCur; NameOfFormulation NumCur;
+ Quantity {
+ { Name u_cur; Value { Local { [ {u_cur} ]; In Region[{BndExt}]; Jacobian JSur; }}}
+ }
+ }
+ { Name NumSol; NameOfFormulation NumSol;
+ Quantity {
+ { Name u_ref; Value { Local { [ f_ref[] ]; In Dom; Jacobian JVol; }}}
+ { Name u_num~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ {u_num} ]; In Region[{Dom,DomPml}]; Jacobian JVol; }}}
+ { Name u_err~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ f_ref[]-{u_num} ]; In Dom; 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; }}}
+ }
+ }
+}
+
+PostOperation{
+ { Name NumNormal; NameOfPostProcessing NumNormal;
+ Operation {
+ Print [u_normal, OnElementsOf Region[{BndExt}], StoreInField (1001), Depth 10];
+// Print [u_normal, OnElementsOf Region[{BndExt}], File "out/u_normal.pos"];
+ }
+ }
+ { Name NumCur; NameOfPostProcessing NumCur;
+ Operation {
+ Print [u_cur, OnElementsOf Region[{BndExt}], StoreInField (1002), Depth 10];
+// Print [u_cur, OnElementsOf Region[{BndExt}], File "out/u_cur.pos"];
+ }
+ }
+ { Name NumSol; NameOfPostProcessing NumSol;
+ Operation {
+ tmp1 = Sprintf("out/solNum_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ tmp2 = Sprintf("out/solErr_%g_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ Print [u_ref, OnElementsOf Dom, File "out/solRef.pos"];
+ Print [u_num~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Region[{Dom,DomPml}], File tmp1];
+ Print [u_err~{BND_TYPE}~{CRN_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Dom, File tmp2];
+ }
+ }
+ { Name Errors; NameOfPostProcessing Errors;
+ Operation {
+ tmp3 = Sprintf("out/errorAbs_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ tmp4 = Sprintf("out/errorRel_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_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_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_TYPE, nPade, thetaPadeInput);
+ tmp6 = Sprintf("out/errorProjRel_%g_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, CRN_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];
+ }
+ }
+}
+
+DefineConstant[
+ R_ = {"NumSol", Name "GetDP/1ResolutionChoices", Visible 1, Choices {"NumNormal", "NumCur", "NumSol"} },
+ P_ = {"NumSol, Errors", Name "GetDP/2PostOperationChoices", Visible 1, Choices {"NumNormal", "NumCur", "NumSol", "Errors", "ProjError"}}
+];