Starting RCot TC (looks good for 2 equal subdomains)

GetDDM/Helmholtz.pro

@@ -124,7 +124,7 @@ FunctionSpace {
           }
         }
       }
-      If (TC_TYPE == 2)
+      If (TC_TYPE == 2 || TC_TYPE == 4)
         For h In {1:NP_OSRC}
           { Name Hgrad_phi~{h}~{i}~{j}; Type Form0;
             BasisFunction {
@@ -148,7 +148,7 @@ Formulation {
         { Name u~{i}; Type Local; NameOfSpace Hgrad_u~{i}; }
         For jj In {0: #myD~{i}()-1}
           j = myD~{i}(jj);
-          If(TC_TYPE == 2)
+          If(TC_TYPE == 2 || TC_TYPE == 4)
             For h In{1:NP_OSRC}
               { Name phi~{h}~{i}~{j}; Type Local;
                 NameOfSpace Hgrad_phi~{h}~{i}~{j}; }
@@ -226,6 +226,28 @@ Formulation {
           EndFor
         EndIf
 
+        If(TC_TYPE == 4) // RCot
+          For jj In {0: #myD~{i}()-1}
+            j = myD~{i}(jj);
+            Galerkin { [ 1/Ld~{j} * Dof{u~{i}}, {u~{i}} ];
+              In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+
+            For h In{1:NP_OSRC}
+              Galerkin { [ +2/Ld~{j} * Dof{phi~{h}~{i}~{j}}, {u~{i}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+              Galerkin { [ -2/Ld~{j} * (1/k[])^2 * Dof{d phi~{h}~{i}~{j}}, {d u~{i}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+
+              Galerkin { [ -1 * Dof{u~{i}}, {phi~{h}~{i}~{j}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+              Galerkin { [ (1-((h*Pi)/(k[]*Ld~{j}))^2) * Dof{phi~{h}~{i}~{j}}, {phi~{h}~{i}~{j}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+              Galerkin { [ -(1/k[])^2 * Dof{d phi~{h}~{i}~{j}}, {d phi~{h}~{i}~{j}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+            EndFor
+          EndFor
+        EndIf
+
         // Bayliss-Turkel absorbing boundary condition
         Galerkin { [ - I[] * kInf[] * Dof{u~{i}} , {u~{i}} ];
           In GammaInf~{i}; Jacobian JSur; Integration I1; }
@@ -246,7 +268,7 @@ Formulation {
           { Name u~{i}; Type Local; NameOfSpace Hgrad_u~{i}; }
           { Name g_out~{i}~{j}; Type Local;
             NameOfSpace Hgrad_g_out~{i}~{j}; }
-          If(TC_TYPE == 2)
+          If(TC_TYPE == 2 || TC_TYPE == 4)
             For h In{1:NP_OSRC}
               { Name phi~{h}~{i}~{j}; Type Local;
                 NameOfSpace Hgrad_phi~{h}~{i}~{j}; }
@@ -290,6 +312,18 @@ Formulation {
             Galerkin { [ 2 * I[] * kInf[] * {u~{i}}, {g_out~{i}~{j}}];
               In TrBndPmlSigma~{i}~{j}; Jacobian JSur; Integration I1;}
           EndIf
+
+          If(TC_TYPE == 4) // RCot
+            Galerkin { [ -2 * 1/Ld~{j} * {u~{i}}, {g_out~{i}~{j}} ];
+              In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+            For h In{1:NP_OSRC}
+              Galerkin { [ -2 * 2/Ld~{j} * {u~{i}}, {g_out~{i}~{j}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+              Galerkin { [ -2 * 2/Ld~{j} * ((h*Pi)/(k[]*Ld~{j}))^2 * {phi~{h}~{i}~{j}}, {g_out~{i}~{j}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+            EndFor
+          EndIf
+
         }
       }
 

GetDDM/cavity2d.pro

+Include "cavity2d_data.geo";
+
+DefineConstant[
+  TC_TYPE = {2, Name "Input/01Transmission condition",
+    Choices {0="Order 0", 1="Order 2", 2="Pade (OSRC)", 4="RCot"}},
+  NP_OSRC = {4, Name "Input/02Number of auxiliary variables (OSRC or RCot)"},
+  PRECONDITIONER = 0,
+  IS_OPEN = {0, Name "Input/03Is cavity open?", Choices{0, 1}},
+  MODE    = {-1, Name "Input/04Mode"},
+  ListOfCuts() = { {0, N_DOM-1} }
+];
+
+Include "cavity2d_src.pro";
+
+Function{
+  // Imaginary value
+  I[] = Complex[0, 1];
+
+  // Domain length
+  Ld~{0} = LX * t2Dom;
+  Ld~{1} = LX * (1-t2Dom);
+
+  // Wavenumber
+  k   = WAVENUMBER;
+  k[] = k;
+
+  // ABC
+  kInf[]    = k;
+  alphaBT[] = 0;
+  betaBT[]  = 0;
+
+  // 0th order TC
+  kIBC[] = k[];
+
+  // 2nd order TC (Gander 2002, pp. 46-47)
+  xsimin      = 0;
+  xsimax      = Pi / LC;
+  deltak[]    = Pi;
+  alphastar[] = I[] * ((k[]^2 - xsimin^2) * (k[]^2 - (k[]-deltak[])^2))^(1/4);
+  betastar[]  = ((xsimax^2 - k[]^2) * ((k[]+deltak[])^2 - k[]^2))^(1/4);
+  a[] = - (alphastar[] * betastar[] - k[]^2) / (alphastar[] + betastar[]);
+  b[] = - 1 / (alphastar[] + betastar[]);
+
+  // Pade TC
+  kepsI        = 0;
+  keps[]       = k[]*(1+kepsI*I[]);
+  theta_branch = Pi/4;
+
+  // No PMLs
+  D[] = 1;
+  E[] = 1;
+}
+
+Group{
+  D() = {};
+  For idom In {0:N_DOM-1}
+    left  = (idom-1)%N_DOM; // left  neighbour
+    right = (idom+1)%N_DOM; // right neighbour
+
+    D() += idom;
+
+    Omega~{idom}    = Region[(idom + 1)];
+    GammaD~{idom}   = Region[{}];
+    GammaD0~{idom}  = Region[(idom + 2*N_DOM + 2)];
+    GammaInf~{idom} = Region[{}];
+    GammaN~{idom}   = Region[{}];
+
+    If(idom == 0)
+      D~{idom} = {right};
+
+      Sigma~{idom}~{right} = Region[(idom + N_DOM + 2)];
+      Sigma~{idom}         = Region[{Sigma~{idom}~{right}}];
+      GammaD~{idom} += Region[(idom + N_DOM + 1)];
+
+      BndSigma~{idom}~{right}    = Region[{}];
+      BndSigma~{idom}            = Region[{BndSigma~{idom}~{right}}];
+      BndGammaInf~{idom}~{right} = Region[{}];
+      BndGammaInf~{idom}         = Region[{BndGammaInf~{idom}~{right}}];
+
+      Pml~{idom}~{right}    = Region[{}];
+      PmlD0~{idom}~{right}  = Region[{}];
+      PmlInf~{idom}~{right} = Region[{}];
+    EndIf
+    If(idom == N_DOM-1)
+      D~{idom} = {left};
+
+      Sigma~{idom}~{left}  = Region[(idom + N_DOM + 1)];
+      Sigma~{idom}         = Region[{Sigma~{idom}~{left}}];
+
+      If(IS_OPEN)
+        GammaInf~{idom} += Region[(idom + N_DOM + 2)];
+      Else
+        GammaD0~{idom}  += Region[(idom + N_DOM + 2)];
+      EndIf
+      BndSigma~{idom}~{left}    = Region[{}];
+      BndSigma~{idom}           = Region[{BndSigma~{idom}~{left}}];
+      BndGammaInf~{idom}~{left} = Region[{}];
+      BndGammaInf~{idom}        = Region[{BndGammaInf~{idom}~{left}}];
+
+      Pml~{idom}~{left}    = Region[{}];
+      PmlD0~{idom}~{left}  = Region[{}];
+      PmlInf~{idom}~{left} = Region[{}];
+    EndIf
+    If(idom >= 1 && idom < N_DOM-1)
+      D~{idom} = {left, right};
+
+      Sigma~{idom}~{left}  = Region[(idom + N_DOM + 1)];
+      Sigma~{idom}~{right} = Region[(idom + N_DOM + 2)];
+      Sigma~{idom}         = Region[{Sigma~{idom}~{left},
+                                     Sigma~{idom}~{right}}];
+
+      BndSigma~{idom}~{left}     = Region[{}];
+      BndSigma~{idom}~{right}    = Region[{}];
+      BndSigma~{idom}            = Region[{BndSigma~{idom}~{left},
+                                           BndSigma~{idom}~{right}}];
+      BndGammaInf~{idom}~{left}  = Region[{}];
+      BndGammaInf~{idom}~{right} = Region[{}];
+      BndGammaInf~{idom}         = Region[{BndGammaInf~{idom}~{left},
+                                           BndGammaInf~{idom}~{right}}];
+
+      Pml~{idom}~{left}     = Region[{}];
+      Pml~{idom}~{right}    = Region[{}];
+      PmlD0~{idom}~{left}   = Region[{}];
+      PmlD0~{idom}~{right}  = Region[{}];
+      PmlInf~{idom}~{left}  = Region[{}];
+      PmlInf~{idom}~{right} = Region[{}];
+    EndIf
+  EndFor
+}
+
+Include "Decomposition.pro";
+Include "Helmholtz.pro";
+Include "Schwarz.pro";

GetDDM/cavity2d_src.pro

+Function{
+  If(MODE > 0)
+    ky     = MODE * Pi / LY;
+    uinc[] = Complex[Sin[ky * Y[]], 0];
+  Else
+    uinc[] = Complex[Sin[ 1 * Pi/LY * Y[]], 0] +
+             Complex[Sin[ 2 * Pi/LY * Y[]], 0] +
+             Complex[Sin[ 3 * Pi/LY * Y[]], 0] +
+             Complex[Sin[ 4 * Pi/LY * Y[]], 0] +
+             Complex[Sin[ 5 * Pi/LY * Y[]], 0] +
+             Complex[Sin[ 6 * Pi/LY * Y[]], 0] +
+             Complex[Sin[ 7 * Pi/LY * Y[]], 0] +
+             Complex[Sin[ 8 * Pi/LY * Y[]], 0] +
+             Complex[Sin[ 9 * Pi/LY * Y[]], 0] +
+
+             Complex[Sin[10 * Pi/LY * Y[]], 0] +
+             Complex[Sin[11 * Pi/LY * Y[]], 0] +
+             Complex[Sin[12 * Pi/LY * Y[]], 0] +
+             Complex[Sin[13 * Pi/LY * Y[]], 0] +
+             Complex[Sin[14 * Pi/LY * Y[]], 0] +
+             Complex[Sin[15 * Pi/LY * Y[]], 0] +
+             Complex[Sin[16 * Pi/LY * Y[]], 0] +
+             Complex[Sin[17 * Pi/LY * Y[]], 0] +
+             Complex[Sin[18 * Pi/LY * Y[]], 0] +
+             Complex[Sin[19 * Pi/LY * Y[]], 0] +
+
+             Complex[Sin[20 * Pi/LY * Y[]], 0] +
+             Complex[Sin[21 * Pi/LY * Y[]], 0] +
+             Complex[Sin[22 * Pi/LY * Y[]], 0] +
+             Complex[Sin[23 * Pi/LY * Y[]], 0] +
+             Complex[Sin[24 * Pi/LY * Y[]], 0] +
+             Complex[Sin[25 * Pi/LY * Y[]], 0] +
+             Complex[Sin[26 * Pi/LY * Y[]], 0] +
+             Complex[Sin[27 * Pi/LY * Y[]], 0] +
+             Complex[Sin[28 * Pi/LY * Y[]], 0] +
+             Complex[Sin[29 * Pi/LY * Y[]], 0] +
+
+             Complex[Sin[30 * Pi/LY * Y[]], 0] +
+             Complex[Sin[31 * Pi/LY * Y[]], 0] +
+             Complex[Sin[32 * Pi/LY * Y[]], 0] +
+             Complex[Sin[33 * Pi/LY * Y[]], 0] +
+             Complex[Sin[34 * Pi/LY * Y[]], 0] +
+             Complex[Sin[35 * Pi/LY * Y[]], 0] +
+             Complex[Sin[36 * Pi/LY * Y[]], 0] +
+             Complex[Sin[37 * Pi/LY * Y[]], 0] +
+             Complex[Sin[38 * Pi/LY * Y[]], 0] +
+             Complex[Sin[39 * Pi/LY * Y[]], 0] +
+
+             Complex[Sin[40 * Pi/LY * Y[]], 0] +
+             Complex[Sin[41 * Pi/LY * Y[]], 0] +
+             Complex[Sin[42 * Pi/LY * Y[]], 0] +
+             Complex[Sin[43 * Pi/LY * Y[]], 0] +
+             Complex[Sin[44 * Pi/LY * Y[]], 0] +
+             Complex[Sin[45 * Pi/LY * Y[]], 0] +
+             Complex[Sin[46 * Pi/LY * Y[]], 0] +
+             Complex[Sin[47 * Pi/LY * Y[]], 0] +
+             Complex[Sin[48 * Pi/LY * Y[]], 0] +
+             Complex[Sin[49 * Pi/LY * Y[]], 0] +
+
+             Complex[Sin[50 * Pi/LY * Y[]], 0] +
+             Complex[Sin[51 * Pi/LY * Y[]], 0] +
+             Complex[Sin[52 * Pi/LY * Y[]], 0] +
+             Complex[Sin[53 * Pi/LY * Y[]], 0] +
+             Complex[Sin[54 * Pi/LY * Y[]], 0] +
+             Complex[Sin[55 * Pi/LY * Y[]], 0] +
+             Complex[Sin[56 * Pi/LY * Y[]], 0] +
+             Complex[Sin[57 * Pi/LY * Y[]], 0] +
+             Complex[Sin[58 * Pi/LY * Y[]], 0] +
+             Complex[Sin[59 * Pi/LY * Y[]], 0] +
+
+             Complex[Sin[60 * Pi/LY * Y[]], 0] +
+             Complex[Sin[61 * Pi/LY * Y[]], 0] +
+             Complex[Sin[62 * Pi/LY * Y[]], 0] +
+             Complex[Sin[63 * Pi/LY * Y[]], 0] +
+             Complex[Sin[64 * Pi/LY * Y[]], 0] +
+             Complex[Sin[65 * Pi/LY * Y[]], 0] +
+             Complex[Sin[66 * Pi/LY * Y[]], 0] +
+             Complex[Sin[67 * Pi/LY * Y[]], 0] +
+             Complex[Sin[68 * Pi/LY * Y[]], 0] +
+             Complex[Sin[69 * Pi/LY * Y[]], 0] +
+
+             Complex[Sin[70 * Pi/LY * Y[]], 0] +
+             Complex[Sin[71 * Pi/LY * Y[]], 0] +
+             Complex[Sin[72 * Pi/LY * Y[]], 0] +
+             Complex[Sin[73 * Pi/LY * Y[]], 0] +
+             Complex[Sin[74 * Pi/LY * Y[]], 0] +
+             Complex[Sin[75 * Pi/LY * Y[]], 0] +
+             Complex[Sin[76 * Pi/LY * Y[]], 0] +
+             Complex[Sin[77 * Pi/LY * Y[]], 0] +
+             Complex[Sin[78 * Pi/LY * Y[]], 0] +
+             Complex[Sin[79 * Pi/LY * Y[]], 0] +
+
+             Complex[Sin[80 * Pi/LY * Y[]], 0] +
+             Complex[Sin[81 * Pi/LY * Y[]], 0] +
+             Complex[Sin[82 * Pi/LY * Y[]], 0] +
+             Complex[Sin[83 * Pi/LY * Y[]], 0] +
+             Complex[Sin[84 * Pi/LY * Y[]], 0] +
+             Complex[Sin[85 * Pi/LY * Y[]], 0] +
+             Complex[Sin[86 * Pi/LY * Y[]], 0] +
+             Complex[Sin[87 * Pi/LY * Y[]], 0] +
+             Complex[Sin[88 * Pi/LY * Y[]], 0] +
+             Complex[Sin[89 * Pi/LY * Y[]], 0] +
+
+             Complex[Sin[90 * Pi/LY * Y[]], 0] +
+             Complex[Sin[91 * Pi/LY * Y[]], 0] +
+             Complex[Sin[92 * Pi/LY * Y[]], 0] +
+             Complex[Sin[93 * Pi/LY * Y[]], 0] +
+             Complex[Sin[94 * Pi/LY * Y[]], 0] +
+             Complex[Sin[95 * Pi/LY * Y[]], 0] +
+             Complex[Sin[96 * Pi/LY * Y[]], 0] +
+             Complex[Sin[97 * Pi/LY * Y[]], 0] +
+             Complex[Sin[98 * Pi/LY * Y[]], 0] +
+             Complex[Sin[99 * Pi/LY * Y[]], 0];
+  EndIf
+}