RCot for non equal domains by adding psi aux. field + separate phi aux. field...

RCot for non equal domains by adding psi aux. field + separate phi aux. field for RCOT and OSRC + seems to work also with 3 subdomains \o/

GetDDM/Helmholtz.pro

@@ -124,7 +124,7 @@ FunctionSpace {
           }
         }
       }
-      If (TC_TYPE == 2 || TC_TYPE == 4)
+      If (TC_TYPE == 2)
         For h In {1:NP_OSRC}
           { Name Hgrad_phi~{h}~{i}~{j}; Type Form0;
             BasisFunction {
@@ -136,6 +136,26 @@ FunctionSpace {
           }
         EndFor
       EndIf
+      If (TC_TYPE == 4)
+        For h In {1:NP_RCOT}
+          { Name Hgrad_phi~{h}~{i}~{j}; Type Form0;
+            BasisFunction {
+              { Name sn; NameOfCoef un; Function BF_Node;
+                Support Region[ {Sigma~{i}~{j}, BndSigmaInf~{i}~{j}, BndSigmaN~{i}~{j}} ];
+                Entity NodesOf[All, Not {GammaD~{i}, GammaD0~{i}}];
+              }
+            }
+          }
+          { Name Hgrad_psi~{h}~{i}~{j}; Type Form0;
+            BasisFunction {
+              { Name sn; NameOfCoef un; Function BF_Node;
+                Support Region[ {Sigma~{i}~{j}, BndSigmaInf~{i}~{j}, BndSigmaN~{i}~{j}} ];
+                Entity NodesOf[All, Not {GammaD~{i}, GammaD0~{i}}];
+              }
+            }
+          }
+        EndFor
+      EndIf
     EndFor
   EndFor
 }
@@ -148,12 +168,18 @@ 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 || TC_TYPE == 4)
+          If(TC_TYPE == 2)
             For h In{1:NP_OSRC}
               { Name phi~{h}~{i}~{j}; Type Local;
                 NameOfSpace Hgrad_phi~{h}~{i}~{j}; }
             EndFor
           EndIf
+          If(TC_TYPE == 4)
+            For h In{1:NP_RCOT}
+              { Name phi~{h}~{i}~{j}; Type Local;
+                NameOfSpace Hgrad_phi~{h}~{i}~{j}; }
+            EndFor
+          EndIf
         EndFor
       }
       Equation {
@@ -232,7 +258,7 @@ Formulation {
             Galerkin { [ 1/Ld~{j} * Dof{u~{i}}, {u~{i}} ];
               In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
 
-            For h In{1:NP_OSRC}
+            For h In{1:NP_RCOT}
               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}} ];
@@ -268,12 +294,20 @@ 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 || TC_TYPE == 4)
+          If(TC_TYPE == 2)
             For h In{1:NP_OSRC}
               { Name phi~{h}~{i}~{j}; Type Local;
                 NameOfSpace Hgrad_phi~{h}~{i}~{j}; }
             EndFor
           EndIf
+          If(TC_TYPE == 4)
+            For h In{1:NP_RCOT}
+              { Name phi~{h}~{i}~{j}; Type Local;
+                NameOfSpace Hgrad_phi~{h}~{i}~{j}; }
+              { Name psi~{h}~{i}~{j}; Type Local;
+                NameOfSpace Hgrad_psi~{h}~{i}~{j}; }
+            EndFor
+          EndIf
         }
         Equation {
           // reverse sign if d_n computed in Omega
@@ -314,16 +348,36 @@ Formulation {
           EndIf
 
           If(TC_TYPE == 4) // RCot
-            Galerkin { [ -2 * 1/Ld~{j} * {u~{i}}, {g_out~{i}~{j}} ];
+            // Constribution of S~{i}~{j}(u~{i})
+            Galerkin { [ -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}} ];
+            For h In{1:NP_RCOT}
+              Galerkin { [ -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}} ];
+              Galerkin { [ -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
 
+            // Constribution of S~{j}~{i}(u~{i})
+            Galerkin { [ -1/Ld~{i} * {u~{i}}, {g_out~{i}~{j}} ];
+              In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+            For h In{1:NP_RCOT}
+              Galerkin { [ -2/Ld~{i} * Dof{psi~{h}~{i}~{j}}, {g_out~{i}~{j}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+              Galerkin { [ +2/Ld~{i} * (1/k[])^2 * Dof{d psi~{h}~{i}~{j}}, {d g_out~{i}~{j}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+            EndFor
+
+            // Equations for psi auxiliary field
+            For h In{1:NP_RCOT}
+              Galerkin { [ -1 * {u~{i}}, {psi~{h}~{i}~{j}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+              Galerkin { [ (1-((h*Pi)/(k[]*Ld~{i}))^2) * Dof{psi~{h}~{i}~{j}}, {psi~{h}~{i}~{j}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+              Galerkin { [ -(1/k[])^2 * Dof{d psi~{h}~{i}~{j}}, {d psi~{h}~{i}~{j}} ];
+                In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }
+            EndFor
+          EndIf
         }
       }
 

GetDDM/cavity2d.pro

@@ -3,13 +3,16 @@ 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)"},
+  NP_AUX = {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} }
 ];
 
+NP_OSRC = NP_AUX;
+NP_RCOT = NP_AUX;
+
 Include "cavity2d_src.pro";
 
 Function{
@@ -17,8 +20,14 @@ Function{
   I[] = Complex[0, 1];
 
   // Domain length
-  Ld~{0} = LX * t2Dom;
-  Ld~{1} = LX * (1-t2Dom);
+  If(N_DOM == 2)
+    Ld~{0} = LX * t2Dom;
+    Ld~{1} = LX * (1-t2Dom);
+  Else
+    For i In {0:N_DOM-1}
+      Ld~{i} = LX/N_DOM;
+    EndFor
+  EndIf
 
   // Wavenumber
   k   = WAVENUMBER;

GetDDM/cavity2d_data.geo

@@ -13,7 +13,7 @@ DefineConstant[
   LZ = {0.25, Name "Input/3Z dimension"},
 
   N_DOM = {2, Name "Input/4Number of subdomains"},
-  t2Dom = 0.5,
+  t2Dom = 0.33,
 
   QUAD = {0, Name "Input/5Quad mesh?", Choices{0, 1}},