add model HelmholtzDDMwithCrosspoints

HelmholtzDDMwithCrosspoints/README.txt

+A non-overlapping HABC-based domain decomposition method with cross-point treatment for 2D Helmholtz problems.
+
+A. Modave, A. Royer, X. Antoine, X. Geuzaine. An optimized Schwarz domain decomposition method with cross-point treatment for time-harmonic acoustic scattering. https://hal.archives-ouvertes.fr/hal-02432422
+
+Model developed by Axel Modave, based on GetDDM codes.
+
+
+Quick start
+-----------
+
+Open `main.pro' with Gmsh.
\ No newline at end of file

HelmholtzDDMwithCrosspoints/circularDom.dat

+BC_Dir = 1;
+BC_Neu = 2;
+
+DefineConstant[
+  BC_SCATT = {BC_Dir, Highlight "Blue",
+    Choices {BC_Dir = "Sound-soft scat. (Dirichlet BC)",
+             BC_Neu = "Sound-hard scat. (Neumann BC)"},
+    Name "Input/01Type of solution"},
+
+  // Geometry
+  R_INT = {1, Min 0.01, Step 0.01, Name "Input/1Model/40Internal radius"},
+  R_EXT = {4, Min 0.01, Step 0.01, Name "Input/1Model/41External radius"},
+
+  // Subdomains
+  Npie_DOM = {5, Min 3, Max 20, Step 1, Name "Input/1Model/50Number of azimutal divisions"},
+  Nlay_DOM = {3, Min 3, Max 20, Step 1, Name "Input/1Model/51Number of radial divisions"},
+  THETA_INC = {0, Min 0., Max 2*Pi, Step 0.1, Name "Input/1Model/52Rotation of pie"},
+  TWIST   = {0, Min 0, Max 0.3, Step 0.1, Name "Input/1Model/52 TWIST"},
+
+  // Frequency
+  WAVENUMBER = {4*Pi, Min 1, Max 24, Step 1, Name "Input/03Wavenumber"},
+  N_LAMBDA   = {10,   Min 5, Max 30, Step 1, Name "Input/04Points per wavelength"}
+];
+
+LAMBDA = (2*Pi)/WAVENUMBER;
+FREQ = WAVENUMBER/(2*Pi);
+LC = LAMBDA/N_LAMBDA;
+
+N_DOM = Npie_DOM*Nlay_DOM;
+R_SCA = R_INT;
+
+For iDom In {0:N_DOM-1}
+  TAG_DOM~{iDom} = 1000 + iDom;
+For iEdge In {0:3}
+  TAG_BND~{iDom}~{iEdge} = 2000 + 4*iDom + iEdge;
+EndFor
+For iCorner In {0:3}
+  TAG_CRN~{iDom}~{iCorner} = 3000 + 4*iDom + iCorner;
+EndFor
+EndFor

HelmholtzDDMwithCrosspoints/circularDom.geo

+Include "circularDom.dat";
+
+Mesh.CharacteristicLengthMin = LC;
+Mesh.CharacteristicLengthMax = LC;
+
+// =================================================================================================
+// Build POINTS, RADIAL LINES, CIRCULAR LINES and SURFACES
+// =================================================================================================
+
+If(Npie_DOM == 1)
+
+  // Circular lines
+  pCenter = newp;
+  Point(pCenter) = {0,0,0};
+  For iLay In {0:Nlay_DOM}
+    radius = R_INT + iLay * (R_EXT-R_INT)/Nlay_DOM;
+    p~{iLay}~{0} = newp;
+    Point(p~{iLay}~{0}) = { radius, 0, 0};
+    p~{iLay}~{1} = newp;
+    Point(p~{iLay}~{1}) = {-radius, 0, 0};
+    l_curve~{iLay}~{0} = newl;
+    Circle(l_curve~{iLay}~{0}) = {p~{iLay}~{0}, pCenter, p~{iLay}~{1}};
+    l_curve~{iLay}~{1} = newl;
+    Circle(l_curve~{iLay}~{1}) = {p~{iLay}~{1}, pCenter, p~{iLay}~{0}};
+  EndFor
+
+  // Surfaces
+  For iLay In {0:Nlay_DOM-1}
+    ll = newll;
+    Line Loop(ll) = {l_curve~{iLay+1}~{0}, l_curve~{iLay+1}~{1}, -l_curve~{iLay}~{0}, -l_curve~{iLay}~{1} };
+    omega~{iLay}~{0} = news;
+    Plane Surface(omega~{iLay}~{0}) = {ll};
+  EndFor
+
+Else
+
+  // Points
+  
+  radiusTwist~{0} = -TWIST;
+  radiusTwist~{1} =  TWIST;
+  radiusTwist~{2} = -TWIST;
+  radiusTwist~{3} =  TWIST;
+  radiusTwist~{4} = -TWIST;
+  radiusTwist~{5} =  TWIST;
+  radiusTwist~{6} = -TWIST;
+  radiusTwist~{7} =  TWIST;
+  
+  pCenter = newp;
+  Point(pCenter) = {0,0,0};
+  For iPie In {0:Npie_DOM-1}
+  For iLay In {0:Nlay_DOM}
+    theta  = iPie * (2*Pi/Npie_DOM) - THETA_INC + radiusTwist~{iLay};
+    radius = R_INT + iLay * (R_EXT-R_INT)/Nlay_DOM;
+    p~{iLay}~{iPie} = newp;
+    Point(p~{iLay}~{iPie}) = {radius*Cos(theta), radius*Sin(theta), 0};
+  EndFor
+  EndFor
+
+  // Radial lines
+  For iPie In {0:Npie_DOM-1}
+  For iLay In {0:Nlay_DOM-1}
+    l_straight~{iLay}~{iPie} = newl;
+    Line(l_straight~{iLay}~{iPie}) = {p~{iLay}~{iPie}, p~{iLay+1}~{iPie}};
+  EndFor
+  EndFor
+
+  // Circular lines
+  For iPie In {0:Npie_DOM-1}
+  For iLay In {0:Nlay_DOM}
+    iPieNext = (iPie+1)%Npie_DOM;
+    l_curve~{iLay}~{iPie} = newl;
+    Circle(l_curve~{iLay}~{iPie}) = {p~{iLay}~{iPie}, pCenter, p~{iLay}~{iPieNext}};
+  EndFor
+  EndFor
+
+  // Surfaces
+  For iPie In {0:Npie_DOM-1}
+  For iLay In {0:Nlay_DOM-1}
+    iPieNext = (iPie+1)%Npie_DOM;
+    ll = newll;
+    Line Loop(ll) = { l_straight~{iLay}~{iPie},      l_curve~{iLay+1}~{iPie},
+                     -l_straight~{iLay}~{iPieNext}, -l_curve~{iLay  }~{iPie} };
+    omega~{iLay}~{iPie} = news;
+    Plane Surface(omega~{iLay}~{iPie}) = {ll};
+  EndFor
+  EndFor
+
+EndIf
+
+// =================================================================================================
+// Generate MESH for each subdomain and reference domain
+// =================================================================================================
+
+// Mesh
+If(StrCmp(OnelabAction, "check")) // only mesh if not in onelab check mode
+  Mesh 2;
+EndIf
+
+For iPie In {0:Npie_DOM-1}
+For iLay In {0:Nlay_DOM-1}
+  iDom = iLay + Nlay_DOM*iPie;
+  iPieNext = (iPie+1)%Npie_DOM;
+
+  Delete Physicals;
+
+If(Npie_DOM == 1)
+  Physical Line(TAG_BND~{iDom}~{1}) = { l_curve~{iLay+1}~{0}, l_curve~{iLay+1}~{1} };  // Up
+  Physical Line(TAG_BND~{iDom}~{3}) = { l_curve~{iLay  }~{0}, l_curve~{iLay  }~{1} };  // Down
+Else
+  Physical Point(TAG_CRN~{iDom}~{0}) = { p~{iLay}~{iPie} };
+  Physical Point(TAG_CRN~{iDom}~{1}) = { p~{iLay+1}~{iPie} };
+  Physical Point(TAG_CRN~{iDom}~{2}) = { p~{iLay+1}~{iPieNext} };
+  Physical Point(TAG_CRN~{iDom}~{3}) = { p~{iLay}~{iPieNext} };
+  Physical Line(TAG_BND~{iDom}~{0}) = { l_straight~{iLay}~{iPie}     };  // Right
+  Physical Line(TAG_BND~{iDom}~{1}) = { l_curve~{iLay+1}~{iPie}      };  // Up
+  Physical Line(TAG_BND~{iDom}~{2}) = {-l_straight~{iLay}~{iPieNext} };  // Left
+  Physical Line(TAG_BND~{iDom}~{3}) = { l_curve~{iLay  }~{iPie}      };  // Down
+EndIf
+
+  Physical Surface(TAG_DOM~{iDom}) = omega~{iLay}~{iPie};
+
+  Printf("Saving subdomain %g...", iDom);
+  Save StrCat(MSH_NAME, Sprintf("_%g.msh", iDom));
+  Printf("Done.");
+
+EndFor
+EndFor
+
+Delete Physicals;
+
+For iPie In {0:Npie_DOM-1}
+For iLay In {0:Nlay_DOM-1}
+  iDom = iLay + Nlay_DOM*iPie;
+  iPieNext = (iPie+1)%Npie_DOM;
+
+If(Npie_DOM == 1)
+If(iLay == 0)
+  Physical Line(TAG_BND~{iDom}~{3}) = { l_curve~{iLay  }~{0}, l_curve~{iLay  }~{1} };  // Down
+EndIf
+If(iLay == Nlay_DOM-1)
+  Physical Line(TAG_BND~{iDom}~{1}) = { l_curve~{iLay+1}~{0}, l_curve~{iLay+1}~{1} };  // Up
+EndIf
+Else
+If(iLay == 0)
+  Physical Line(TAG_BND~{iDom}~{3}) = { l_curve~{iLay  }~{iPie} };  // Down
+EndIf
+If(iLay == Nlay_DOM-1)
+  Physical Line(TAG_BND~{iDom}~{1}) = { l_curve~{iLay+1}~{iPie} };  // Up
+EndIf
+EndIf
+
+  Physical Surface(TAG_DOM~{iDom}) = omega~{iLay}~{iPie};
+
+EndFor
+EndFor
+
+Printf("Saving domain");
+Save StrCat(MSH_NAME, Sprintf(".msh"));
+Printf("Done.");

HelmholtzDDMwithCrosspoints/circularDom.pro

+Include "circularDom.dat";
+
+// ====================================================================================================
+// FUNCTIONS
+// ====================================================================================================
+
+Function{
+  R[] = Sqrt[X[]*X[]+Y[]*Y[]];
+  f_inc[] = Complex[Cos[WAVENUMBER*X[]], Sin[WAVENUMBER*X[]]];
+  dfdx_inc[] = WAVENUMBER * Complex[-Sin[WAVENUMBER*X[]], Cos[WAVENUMBER*X[]]];
+  
+  // Dirichlet
+If(BC_SCATT == BC_Dir)
+  f_ref[] = AcousticFieldSoftCylinder[XYZ[]]{WAVENUMBER, R_SCA};
+  f_dir[] = -f_inc[];
+EndIf
+
+  // Neumann
+If(BC_SCATT == BC_Neu)
+  f_ref[] = AcousticFieldHardCylinder[XYZ[]]{WAVENUMBER, R_SCA, 0, 0, 50};
+  f_neu[] = dfdx_inc[] * (X[]/R[]);
+EndIf
+}
+
+// ====================================================================================================
+// LISTS ABC/TBC
+// ====================================================================================================
+
+// Is there an ABC for main domain? (1=yes, 0=no)
+isEdgeABC~{0} = 1;
+isEdgeABC~{1} = 1;
+isEdgeABC~{2} = 1;
+isEdgeABC~{3} = 1;
+isEdgeRad~{0} = 1;
+isEdgeRad~{1} = 1;
+isEdgeRad~{2} = 1;
+isEdgeRad~{3} = 1;
+
+ListOfSubdom = {};
+For iPie In {0:Npie_DOM-1}
+For iLay In {0:Nlay_DOM-1}
+  iDom = iLay + Nlay_DOM*iPie;
+
+If (iDom % MPI_Size == MPI_Rank)
+  ListOfSubdom += iDom;
+EndIf
+
+  iPiePrev = (iPie+Npie_DOM-1)%Npie_DOM;
+  iPieNext = (iPie+1)%Npie_DOM;
+
+  iDomNeigh~{iDom}~{0} = iLay + Nlay_DOM*iPiePrev;  // Right
+  iDomNeigh~{iDom}~{1} = (iLay+1) + Nlay_DOM*iPie;  // Top
+  iDomNeigh~{iDom}~{2} = iLay + Nlay_DOM*iPieNext;  // Left
+  iDomNeigh~{iDom}~{3} = (iLay-1) + Nlay_DOM*iPie;  // Down
+  iEdgeNeigh~{iDom}~{0} = 2;
+  iEdgeNeigh~{iDom}~{1} = 3;
+  iEdgeNeigh~{iDom}~{2} = 0;
+  iEdgeNeigh~{iDom}~{3} = 1;
+  isEdgeRad~{iDom}~{0} = 0;
+  isEdgeRad~{iDom}~{1} = 1;
+  isEdgeRad~{iDom}~{2} = 0;
+  isEdgeRad~{iDom}~{3} = 1;
+
+  // Is there an ABC or TBC? (1=yes, 0=no)
+  For iEdge In {0:3}
+    isEdgeABC~{iDom}~{iEdge} = 0;
+    isEdgeTBC~{iDom}~{iEdge} = 1;
+  EndFor
+  If(iLay == 0)
+    isEdgeTBC~{iDom}~{3} = 0;
+  EndIf
+  If(iLay == Nlay_DOM-1)
+    isEdgeABC~{iDom}~{1} = 1;
+    isEdgeTBC~{iDom}~{1} = 0;
+    isEdgeRad~{iDom}~{1} = isEdgeRad~{1};
+  EndIf
+  If(Npie_DOM == 1)
+    isEdgeTBC~{iDom}~{0} = 0;
+    isEdgeTBC~{iDom}~{2} = 0;
+  EndIf
+
+  // List of Edges for Absorbing/Transmission Boundary Condition
+  ListOfEdgesWithABC~{iDom} = {};
+  ListOfEdgesWithTBC~{iDom} = {};
+  For iEdge In {0:3}
+    If(isEdgeABC~{iDom}~{iEdge} == 1)
+      ListOfEdgesWithABC~{iDom} += iEdge;
+    EndIf
+    If(isEdgeTBC~{iDom}~{iEdge} == 1)
+      ListOfEdgesWithTBC~{iDom} += iEdge;
+    EndIf
+  EndFor
+
+EndFor
+EndFor
+
+// ====================================================================================================
+// GROUPS
+// ====================================================================================================
+
+Group {
+
+// === REFERENCE CASE ===
+
+  Omega = Region[{}];
+  GammaD = Region[{}];
+  GammaD0 = Region[{}];
+  GammaN = Region[{}];
+  GammaN0 = Region[{}];
+  GammaPoint = Region[{}];
+  For iEdge In {0:3}
+    SigmaMain~{iEdge} = Region[{}];
+    SigmaMainSides~{iEdge} = Region[{}];
+    SigmaMainSide~{iEdge}~{0} = Region[{}];
+    SigmaMainSide~{iEdge}~{1} = Region[{}];
+  EndFor
+  For iCorner In {0:3}
+    CornerMain~{iCorner} = Region[{}];
+  EndFor
+
+For iPie In {0:Npie_DOM-1}
+For iLay In {0:Nlay_DOM-1}
+  iDom = iLay + Nlay_DOM*iPie;
+
+  // DOMAIN
+  Omega += Region[TAG_DOM~{iDom}];
+
+  // BOUNDARY (Dirichlet/Neumann BC)
+  If(iLay == 0)
+  If(BC_SCATT == BC_Dir)
+    GammaD += Region[TAG_BND~{iDom}~{3}];
+  EndIf
+  If(BC_SCATT == BC_Neu)
+    GammaN += Region[TAG_BND~{iDom}~{3}];
+  EndIf
+  EndIf
+
+  // BOUNDARY (Absorbing BC)
+  If(iLay == Nlay_DOM-1)
+    SigmaMain~{1} += Region[{TAG_BND~{iDom}~{1}}];
+  EndIf
+
+EndFor
+EndFor
+
+  SigmaMain = Region[{ SigmaMain~{0}, SigmaMain~{1}, SigmaMain~{2}, SigmaMain~{3} }];
+  CornerMain = Region[{ CornerMain~{0}, CornerMain~{1}, CornerMain~{2}, CornerMain~{3} }];
+
+// === DDM CASE ===
+
+For iPie In {0:Npie_DOM-1}
+For iLay In {0:Nlay_DOM-1}
+  iDom = iLay + Nlay_DOM*iPie;
+
+  // DOMAIN
+  Omega~{iDom} = Region[TAG_DOM~{iDom}];
+
+  // SOURCE
+  GammaPoint~{iDom} = Region[{}];
+  
+  // BOUNDARY (Dirichlet/Neumann BC)
+  GammaD0~{iDom} = Region[{}];
+  GammaN0~{iDom} = Region[{}];
+  GammaD~{iDom} = Region[{}];
+  GammaN~{iDom} = Region[{}];
+  If(iLay == 0)
+  If(BC_SCATT == BC_Dir)
+    GammaD~{iDom} = Region[ TAG_BND~{iDom}~{3} ];  // Down
+  EndIf
+  If(BC_SCATT == BC_Neu)
+    GammaN~{iDom} = Region[ TAG_BND~{iDom}~{3} ];  // Down
+  EndIf
+  EndIf
+
+  // BOUNDARY (Absorbing BC + Transmission BC)
+  For iEdge In {0:3}
+    Sigma~{iDom}~{iEdge}         = Region[{}];
+    SigmaSide~{iDom}~{iEdge}~{0} = Region[{}];
+    SigmaSide~{iDom}~{iEdge}~{1} = Region[{}];
+    GammaD~{iDom}~{iEdge} = Region[{}];
+    GammaD0~{iDom}~{iEdge} = Region[{}];
+    GammaN~{iDom}~{iEdge} = Region[{}];
+    GammaN0~{iDom}~{iEdge} = Region[{}];
+  If((isEdgeABC~{iDom}~{iEdge} == 1) || (isEdgeTBC~{iDom}~{iEdge} == 1))
+    iCornerSide~{0} = iEdge;
+    iCornerSide~{1} = (iEdge+1)%4;
+    iEdgeSide~{0} = (iEdge+3)%4;
+    iEdgeSide~{1} = (iEdge+1)%4;
+    
+    Sigma~{iDom}~{iEdge}         = Region[{TAG_BND~{iDom}~{iEdge}}];
+    SigmaSide~{iDom}~{iEdge}~{0} = Region[{TAG_CRN~{iDom}~{iCornerSide~{0}}}];
+    SigmaSide~{iDom}~{iEdge}~{1} = Region[{TAG_CRN~{iDom}~{iCornerSide~{1}}}];
+    
+  For iSide In {0:1}
+  If((isEdgeABC~{iDom}~{iEdgeSide~{iSide}} == 0) && (isEdgeTBC~{iDom}~{iEdgeSide~{iSide}} == 0))
+  If(BC_SCATT == BC_Dir)
+    GammaD~{iDom}~{iEdge} = Region[{TAG_CRN~{iDom}~{iCornerSide~{iSide}}}];
+  EndIf
+  If(BC_SCATT == BC_Neu)
+    GammaN~{iDom}~{iEdge} = Region[{TAG_CRN~{iDom}~{iCornerSide~{iSide}}}];
+  EndIf
+  EndIf
+  EndFor
+  
+  EndIf
+  EndFor
+  Sigma~{iDom} = Region[{ Sigma~{iDom}~{0}, Sigma~{iDom}~{1}, Sigma~{iDom}~{2}, Sigma~{iDom}~{3} }];
+
+  // CORNER (Absorbing BC + Transmission BC)
+  For iCorner In {0:3}
+    Corner~{iDom}~{iCorner} = Region[{TAG_CRN~{iDom}~{iCorner}}];
+  EndFor
+  Corner~{iDom} = Region[{ Corner~{iDom}~{0}, Corner~{iDom}~{1}, Corner~{iDom}~{2}, Corner~{iDom}~{3} }];
+EndFor
+EndFor
+
+}

HelmholtzDDMwithCrosspoints/ddmHelmholtzCrossPoints.pro

+CROSSPOINT_NO   = 0;
+CROSSPOINT_Full = 2;
+
+DefineConstant[
+  nPadeABC = {0, Choices {0, 1, 2, 3, 4, 5, 6}, Highlight "Yellow", // 6
+    Name "Input/3ABC and TBC/3Pade ABC: Number of fields"},
+  thetaPadeInputABC = {0, Min 0, Max 0.5, Step 0.05, // 0.3
+    Name "Input/3ABC and TBC/4Pade ABC: Rotation of branch cut"},
+  nPadeTBC = {3, Choices {0, 2, 4, 6}, Highlight "Yellow",// Loop 1,
+    Name "Input/3ABC and TBC/5Pade TBC: Number of fields"},
+  thetaPadeInputTBC = {0.3, Min 0.1, Max 0.5, Step 0.1,
+    Name "Input/3ABC and TBC/6Pade TBC: Rotation of branch cut"},
+  stretchPadeInputTBC = {1, Choices {0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50},
+    Name "Input/3ABC and TBC/7Pade ABC: Stretch of branch cut"},
+  CROSSPOINT_TREAT = {CROSSPOINT_Full, Highlight "Blue", //Loop 2,
+    Choices {CROSSPOINT_NO   = "Nothing",
+             CROSSPOINT_Full = "Compatibility everywhere"},
+    Name "Input/3ABC and TBC/8Compatibility Cross-Points"}
+];
+
+// ====================================================================================================
+// FUNCTIONS
+// ====================================================================================================
+
+DefineConstant[ R_SCA, f_sou ];
+For iDom In {0:N_DOM-1}
+  DefineConstant[ f_sou~{iDom} ];
+EndFor
+
+Function{
+
+  DefineFunction[ f_dir, f_neu ];
+
+// === Imaginary unit and wavenumber
+
+  I[] = Complex[0, 1];
+If (FLAG_PBM == PBM_MARMOUSI)
+  k[] = 2*Pi*FREQ/cData[];
+Else
+  k[] = 2*Pi*FREQ;
+EndIf
+
+// === Functions for the Pade ABC and Pade TBC
+
+  mPadeABC = 2*nPadeABC+1;
+  mPadeTBC = 2*nPadeTBC+1;
+
+  For i In {1:nPadeABC}
+    cPadeABC~{i} = Tan[i*Pi/mPadeABC]^2;
+  EndFor
+  For i In {1:nPadeTBC}
+    cPadeTBC~{i} = Tan[i*Pi/mPadeTBC]^2;
+  EndFor
+
+  For iEdge In {0:3}
+    kEps~{iEdge}[] = k[];
+  If(isEdgeRad~{iEdge} == 1)
+    //alphaBT[] = 0.5/R_EXT - (1/R_EXT^2) / (8*(1/R_EXT - I[]*k));
+    //betaBT[] = 0.5/(1/R_EXT - I[]*k);
+
+    paramCurv1~{iEdge}[] = - 0.5 * (1/R[]) + (1/R[])^2 / (8*((1/R[]) - I[]*k[]));
+    paramCurv2~{iEdge}[] = - 0.5 / ((1/R[]) - I[]*k[]);
+    //paramCurv2~{iEdge}[] = - 0.5 / (R[]*k[]*k[]);
+  Else
+    paramCurv1~{iEdge}[] = 0;
+    paramCurv2~{iEdge}[] = 0;
+  EndIf
+  EndFor
+
+  // Gander/Magoules/Nataf OO0
+  kMin[] = Pi/7.5;
+  kMinus[] = k[] - kMin[];
+  kPlus[]  = k[] + kMin[];
+  kMax[] = Pi/LC;
+  //alphaPadeTBC[] = Sqrt[Sqrt[(kMax[]*kMax[]-k[]*k[]) * (k[]*k[]-kMinus[]*kMinus[])]]/k[] * Complex[Cos[Pi/4],Sin[Pi/4]];
+
+  alphaPadeABC[] = Complex[Cos[thetaPadeInputABC*Pi/2.],Sin[thetaPadeInputABC*Pi/2.]];
+  alphaPadeTBC[] = Complex[stretchPadeInputTBC*Cos[thetaPadeInputTBC*Pi/2.],stretchPadeInputTBC*Sin[thetaPadeInputTBC*Pi/2.]];
+
+  For iDom In {0:N_DOM-1}
+  For iEdge In {0:3}
+  If((isEdgeABC~{iDom}~{iEdge} == 1) && (isEdgeRad~{iDom}~{iEdge} == 1))
+    paramCurv1~{iDom}~{iEdge}[] = - 0.5 * (1/R[]) + (1/R[])^2 / (8*((1/R[]) - I[]*k[]));
+    paramCurv2~{iDom}~{iEdge}[] = - 0.5 / ((1/R[]) - I[]*k[]);
+    //paramCurv2~{iDom}~{iEdge}[] = - 0.5/(R[]*k[]*k[]);
+  Else
+    paramCurv1~{iDom}~{iEdge}[] = 0;
+    paramCurv2~{iDom}~{iEdge}[] = 0;
+  EndIf
+  If((isEdgeABC~{iDom}~{iEdge} == 0) && (isEdgeRad~{iDom}~{iEdge} == 1))
+    kEps~{iDom}~{iEdge}[] = k[] + I[] * 0.4 * k[]^(1/3) * (1/R[])^(2/3);  // or 0.6
+  //ElseIf((isEdgeABC~{iDom}~{iEdge} == 0) && (isEdgeRad~{iDom}~{iEdge} == 0))
+  //  kEps~{iDom}~{iEdge}[] = k[] + I[] * k[]/4;
+  Else
+    kEps~{iDom}~{iEdge}[] = k[];
+  EndIf
+  If(isEdgeTBC~{iDom}~{iEdge} == 1)
+    nPade~{iDom}~{iEdge} = nPadeTBC;
+    mPade~{iDom}~{iEdge} = mPadeTBC;
+    For i In {1:nPadeTBC}
+      cPade~{i}~{iDom}~{iEdge} = cPadeTBC~{i};
+    EndFor
+    alphaPade~{iDom}~{iEdge}[] = alphaPadeTBC[];
+  Else
+    nPade~{iDom}~{iEdge} = nPadeABC;
+    mPade~{iDom}~{iEdge} = mPadeABC;
+    For i In {1:nPadeABC}
+      cPade~{i}~{iDom}~{iEdge} = cPadeABC~{i};
+    EndFor
+    alphaPade~{iDom}~{iEdge}[] = alphaPadeABC[];
+  EndIf
+  EndFor
+  EndFor
+
+// === Numerotation of fields
+
+  index = 1;
+  For iDom In {0:N_DOM-1}
+    SaveSOL~{iDom} = index; index = index+1;
+  For iEdge In {0:3}
+    FIELD~{iDom}~{iEdge} = index; index = index+1;
+
+  iEdgeSide~{0} = (iEdge+3)%4;
+  iEdgeSide~{1} = (iEdge+1)%4;
+  For iSide In {0:1}
+  For i In {1:nPade~{iDom}~{iEdgeSide~{iSide}}}
+    FIELD~{i}~{iDom}~{iEdge}~{iSide} = index; index = index+1;
+  EndFor
+  EndFor
+  EndFor
+  EndFor
+
+// === Connectivity of fields
+
+  ListOfFields          = {}; // Tableau flottant contenant les fields
+  ListOfConnectedFields = {}; // fields connected to my fields
+
+  For i0 In {0:#ListOfSubdom()-1}
+  iDom  = ListOfSubdom(i0);
+  For i1 In {0:#ListOfEdgesWithTBC~{iDom}()-1}
+  iEdge = ListOfEdgesWithTBC~{iDom}(i1);
+  iEdgeSide~{0} = (iEdge+3)%4;
+  iEdgeSide~{1} = (iEdge+1)%4;
+  iDomNeigh  = iDomNeigh~{iDom}~{iEdge};
+  iEdgeNeigh = iEdgeNeigh~{iDom}~{iEdge};
+  iSideNeigh~{0} = 1;
+  iSideNeigh~{1} = 0;
+
+    g_in~{iDom}~{iEdge}[]  = ComplexScalarField[XYZ[]]{ FIELD~{iDomNeigh}~{iEdgeNeigh}() };
+    ListOfFields          += FIELD~{iDom}~{iEdge};
+    ListOfConnectedFields += { 1, FIELD~{iDomNeigh}~{iEdgeNeigh}() };
+
+  For iSide In {0:1}
+  If (CROSSPOINT_TREAT == CROSSPOINT_Full)
+  For j In {1:nPade~{iDom}~{iEdgeSide~{iSide}}}
+    g_in~{j}~{iDom}~{iEdge}~{iSide}[]  = ComplexScalarField[XYZ[]]{ FIELD~{j}~{iDomNeigh}~{iEdgeNeigh}~{iSideNeigh~{iSide}}() };
+    ListOfFields                      += FIELD~{j}~{iDom}~{iEdge}~{iSide};
+    ListOfConnectedFields             += { 1, FIELD~{j}~{iDomNeigh}~{iEdgeNeigh}~{iSideNeigh~{iSide}}() };
+  EndFor
+  EndIf
+  EndFor
+
+  EndFor
+  EndFor
+
+// === Reference solution
+
+  FIELD_REF = index+1; index = index+1;
+  If(FLAG_ERROR == ERROR_ANA)
+    f_refErr[] = f_ref[];
+  EndIf
+  If(FLAG_ERROR == ERROR_NUM)
+    f_refErr[] = ComplexScalarField[XYZ[]]{ FIELD_REF() };
+  EndIf
+
+}
+
+// ====================================================================================================
+// CONSTRAINTS
+// ====================================================================================================
+
+Constraint{
+  { Name Dirichlet_u;
+    Case {
+      { Region GammaD;  Value f_dir[]; }
+      { Region GammaD0; Value 0; }
+    }
+  }
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+  { Name Dirichlet_u~{iDom};
+    Case {
+      { Region GammaD~{iDom};  Value $PhysicalSource ? f_dir[] : 0.; }
+      { Region GammaD0~{iDom}; Value 0.; }
+    }
+  }
+For iEdge In {0:3}
+For i In {1:nPade~{iDom}~{iEdge}}
+  { Name Dirichlet_u~{i}~{iDom}~{iEdge};
+    Case {
+      { Region GammaD~{iDom}~{iEdge};  Value $PhysicalSource ? f_dir[] : 0.; }   // SigmaSide~{iDom}~{iEdge}~{0}, SigmaSide~{iDom}~{iEdge}~{1}
+      { Region GammaD0~{iDom}~{iEdge}; Value 0.; }
+    }
+  }
+EndFor
+EndFor
+EndFor
+}
+
+// ====================================================================================================
+// FUNCTION SPACE
+// ====================================================================================================
+
+FunctionSpace {
+
+  { Name Href; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {Omega, SigmaMain, CornerMain, GammaPoint, GammaN} ]; Entity NodesOf[All]; }}
+    Constraint {{ NameOfCoef un; EntityType NodesOf; NameOfConstraint Dirichlet_u; }}
+  }
+For iEdge In {0:3}
+  { Name Href_du~{iEdge}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {SigmaMain~{iEdge}, SigmaMainSide~{iEdge}~{0}, SigmaMainSide~{iEdge}~{1}} ]; Entity NodesOf[All]; }}
+  }
+For i In {1:nPadeABC}
+  { Name Href~{i}~{iEdge}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {SigmaMain~{iEdge}, SigmaMainSide~{iEdge}~{0}, SigmaMainSide~{iEdge}~{1}} ]; Entity NodesOf[All]; }}
+  }
+For iSide In {0:1}
+  { Name Href_du~{i}~{iEdge}~{iSide}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {SigmaMainSide~{iEdge}~{iSide}} ]; Entity NodesOf[All]; }}
+  }
+EndFor
+EndFor
+EndFor
+For iCorner In {0:3}
+For i In {1:nPadeABC}
+For j In {1:nPadeABC}
+  { Name Href~{i}~{j}~{iCorner}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {CornerMain~{iCorner}} ]; Entity NodesOf[All]; }}
+  }
+EndFor
+EndFor
+EndFor
+
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+
+  { Name Hddm_u~{iDom}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {Omega~{iDom}, Sigma~{iDom}, Corner~{iDom}, GammaPoint~{iDom}, GammaN~{iDom}} ]; Entity NodesOf[All]; }}
+    Constraint {{ NameOfCoef un; EntityType NodesOf; NameOfConstraint Dirichlet_u~{iDom}; }}
+  }
+For iEdge In {0:3}
+  { Name Hddm_du~{iDom}~{iEdge}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {Sigma~{iDom}~{iEdge}, SigmaSide~{iDom}~{iEdge}~{0}, SigmaSide~{iDom}~{iEdge}~{1}} ]; Entity NodesOf[All]; }} // , Not {GammaD~{iDom}}
+  }
+For i In {1:nPade~{iDom}~{iEdge}}
+  { Name Hddm_u~{i}~{iDom}~{iEdge}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {Sigma~{iDom}~{iEdge}, SigmaSide~{iDom}~{iEdge}~{0}, SigmaSide~{iDom}~{iEdge}~{1}} ]; Entity NodesOf[All]; }}
+//    Constraint {{ NameOfCoef un; EntityType NodesOf; NameOfConstraint Dirichlet_u~{i}~{iDom}~{iEdge}; }}
+  }
+For iSide In {0:1}
+  { Name Hddm_du~{i}~{iDom}~{iEdge}~{iSide}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {SigmaSide~{iDom}~{iEdge}~{iSide}} ]; Entity NodesOf[All]; }}
+  }
+EndFor
+EndFor
+EndFor
+For iCorner In {0:3}
+iEdgeSide~{0} = (iCorner+3)%4;
+iEdgeSide~{1} = iCorner;
+For i In {1:nPade~{iDom}~{iEdgeSide~{0}}}
+For j In {1:nPade~{iDom}~{iEdgeSide~{1}}}
+  { Name Hddm_u~{i}~{j}~{iDom}~{iCorner}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {Corner~{iDom}~{iCorner}} ]; Entity NodesOf[All]; }}
+  }
+EndFor
+EndFor
+EndFor
+
+For iEdge In {0:3}
+iEdgeSide~{0} = (iEdge+3)%4;
+iEdgeSide~{1} = (iEdge+1)%4;
+  { Name Hddm_g_out~{iDom}~{iEdge}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {Sigma~{iDom}~{iEdge}, SigmaSide~{iDom}~{iEdge}~{0}, SigmaSide~{iDom}~{iEdge}~{1}} ]; Entity NodesOf[All]; }} // , Not {GammaD~{iDom}}
+  }
+For iSide In {0:1}
+For i In {1:nPade~{iDom}~{iEdgeSide~{iSide}}}
+  { Name Hddm_g_out~{i}~{iDom}~{iEdge}~{iSide}; Type Form0;
+    BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[ {SigmaSide~{iDom}~{iEdge}~{iSide}} ]; Entity NodesOf[All]; }}
+  }
+EndFor
+EndFor
+EndFor
+
+EndFor
+}
+
+
+// ====================================================================================================
+// FORMULATION
+// ====================================================================================================
+
+Formulation {
+
+// ====================================================================================================
+// Reference SOLVER without DDM
+// ====================================================================================================
+
+  { Name FormRef; Type FemEquation;
+    Quantity {
+      { Name uRef; Type Local; NameOfSpace Href; }
+    For iEdge In {0:3}
+      { Name duRef~{iEdge}; Type Local; NameOfSpace Href_du~{iEdge}; }
+    For i In {1:nPadeABC}
+      { Name uRef~{i}~{iEdge}; Type Local; NameOfSpace Href~{i}~{iEdge}; }
+    For iSide In {0:1}
+      { Name duRef~{i}~{iEdge}~{iSide}; Type Local; NameOfSpace Href_du~{i}~{iEdge}~{iSide}; }
+    EndFor
+    EndFor
+    EndFor
+    For iCorner In {0:3}
+    For i In {1:nPadeABC}
+    For j In {1:nPadeABC}
+      { Name uRef~{i}~{j}~{iCorner}; Type Local; NameOfSpace Href~{i}~{j}~{iCorner}; }
+    EndFor
+    EndFor
+    EndFor
+    }
+    Equation {
+
+// Field on SURFACE
+
+      Galerkin { [ Dof{d uRef}        , {d uRef} ]; In Omega;             Jacobian JSur; Integration I1; }
+      Galerkin { [ - k[]^2 * Dof{uRef}  , {uRef} ]; In Omega;             Jacobian JSur; Integration I1; }
+      Galerkin { [ - f_sou              , {uRef} ]; In GammaPoint;        Jacobian JSur; Integration I1; }
+      Galerkin { [ - f_neu[]            , {uRef} ]; In GammaN;            Jacobian JLin; Integration I1; }
+    For iEdge In {0:3}
+      Galerkin { [ - Dof{duRef~{iEdge}} , {uRef} ]; In SigmaMain~{iEdge}; Jacobian JLin; Integration I1; }
+    EndFor
+
+// DtN on EDGES
+
+    For iEdge In {0:3}
+      Galerkin { [   Dof{duRef~{iEdge}}                                                          , {duRef~{iEdge}} ]; In SigmaMain~{iEdge}; Jacobian JLin; Integration I1; }
+      Galerkin { [ - paramCurv1~{iEdge}[] * Dof{uRef}                                            , {duRef~{iEdge}} ]; In SigmaMain~{iEdge}; Jacobian JLin; Integration I1; }
+      Galerkin { [ - paramCurv2~{iEdge}[] * Dof{d uRef}                                        , {d duRef~{iEdge}} ]; In SigmaMain~{iEdge}; Jacobian JLin; Integration I1; }
+      Galerkin { [ - I[]*k[]*alphaPadeABC[] * Dof{uRef}                                          , {duRef~{iEdge}} ]; In SigmaMain~{iEdge}; Jacobian JLin; Integration I1; }
+    For i In {1:nPadeABC}
+      Galerkin { [ - I[]*k[]*alphaPadeABC[] * 2./mPadeABC * cPadeABC~{i} * Dof{uRef}             , {duRef~{iEdge}} ]; In SigmaMain~{iEdge}; Jacobian JLin; Integration I1; }
+      Galerkin { [ - I[]*k[]*alphaPadeABC[] * 2./mPadeABC * cPadeABC~{i} * Dof{uRef~{i}~{iEdge}} , {duRef~{iEdge}} ]; In SigmaMain~{iEdge}; Jacobian JLin; Integration I1; }
+    EndFor
+    EndFor
+
+// Fields on EDGES
+
+    For iEdge In {0:3}
+    iCornerSide~{0} = iEdge;
+    iCornerSide~{1} = (iEdge+1)%4;
+    For i In {1:nPadeABC}
+
+      Galerkin { [ Dof{d uRef~{i}~{iEdge}}                                         , {d uRef~{i}~{iEdge}} ]; In SigmaMain~{iEdge};                Jacobian JLin; Integration I1; }
+      Galerkin { [ - kEps~{iEdge}[]^2*(alphaPadeABC[]^2*cPadeABC~{i}+1) * Dof{uRef~{i}~{iEdge}} , {uRef~{i}~{iEdge}} ]; In SigmaMain~{iEdge};                Jacobian JLin; Integration I1; }
+      Galerkin { [ - kEps~{iEdge}[]^2*alphaPadeABC[]^2*(cPadeABC~{i}+1) * Dof{uRef}             , {uRef~{i}~{iEdge}} ]; In SigmaMain~{iEdge};                Jacobian JLin; Integration I1; }
+    For iSide In {0:1}
+      Galerkin { [ - Dof{duRef~{i}~{iEdge}~{iSide}}                                  , {uRef~{i}~{iEdge}} ]; In CornerMain~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    EndFor
+
+    EndFor
+    EndFor
+
+// DtN on EDGES
+
+    For iEdge In {0:3}
+    iCornerSide~{0} = iEdge;
+    iCornerSide~{1} = (iEdge+1)%4;
+    For i In {1:nPadeABC}
+
+    For iSide In {0:1}
+      Galerkin { [   Dof{duRef~{i}~{iEdge}~{iSide}}                                                                , {duRef~{i}~{iEdge}~{iSide}} ]; In CornerMain~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+      Galerkin { [ - I[]*k[]*alphaPadeABC[] * Dof{uRef~{i}~{iEdge}}                                                , {duRef~{i}~{iEdge}~{iSide}} ]; In CornerMain~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    For j In {1:nPadeABC}
+      Galerkin { [ - I[]*k[]*alphaPadeABC[] * 2./mPadeABC * cPadeABC~{j} * Dof{uRef~{i}~{iEdge}}                   , {duRef~{i}~{iEdge}~{iSide}} ]; In CornerMain~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    If (iSide == 0)
+      Galerkin { [ - I[]*k[]*alphaPadeABC[] * 2./mPadeABC * cPadeABC~{j} * Dof{uRef~{j}~{i}~{iCornerSide~{iSide}}} , {duRef~{i}~{iEdge}~{iSide}} ]; In CornerMain~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    Else
+      Galerkin { [ - I[]*k[]*alphaPadeABC[] * 2./mPadeABC * cPadeABC~{j} * Dof{uRef~{i}~{j}~{iCornerSide~{iSide}}} , {duRef~{i}~{iEdge}~{iSide}} ]; In CornerMain~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    EndIf
+    EndFor
+    EndFor
+
+    EndFor
+    EndFor
+
+// Fields on CORNERS
+
+    For iCorner In {0:3}
+    iEdgeSide~{0} = (iCorner+3)%4;
+    iEdgeSide~{1} = iCorner;
+    For i In {1:nPadeABC}
+    For j In {1:nPadeABC}
+      Galerkin { [ Dof{uRef~{i}~{j}~{iCorner}}                                                                                                             , {uRef~{i}~{j}~{iCorner}} ]; In CornerMain~{iCorner}; Jacobian JVol; Integration I1; }
+      Galerkin { [ alphaPadeABC[]^2*(cPadeABC~{j}+1) * Dof{uRef~{i}~{iEdgeSide~{0}}} / (alphaPadeABC[]^2*cPadeABC~{i} + alphaPadeABC[]^2*cPadeABC~{j} + 1) , {uRef~{i}~{j}~{iCorner}} ]; In CornerMain~{iCorner}; Jacobian JVol; Integration I1; }
+      Galerkin { [ alphaPadeABC[]^2*(cPadeABC~{i}+1) * Dof{uRef~{j}~{iEdgeSide~{1}}} / (alphaPadeABC[]^2*cPadeABC~{i} + alphaPadeABC[]^2*cPadeABC~{j} + 1) , {uRef~{i}~{j}~{iCorner}} ]; In CornerMain~{iCorner}; Jacobian JVol; Integration I1; }
+    EndFor
+    EndFor
+    EndFor
+
+    }
+  }
+
+// ====================================================================================================
+// SOLVER for subdomain 'iDom'
+// ====================================================================================================
+
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+
+  { Name FormVol~{iDom}; Type FemEquation;
+    Quantity {
+      { Name u~{iDom}; Type Local; NameOfSpace Hddm_u~{iDom}; }
+    For iEdge In {0:3}
+      { Name du~{iDom}~{iEdge}; Type Local; NameOfSpace Hddm_du~{iDom}~{iEdge}; }
+    For i In {1:nPade~{iDom}~{iEdge}}
+      { Name u~{i}~{iDom}~{iEdge}; Type Local; NameOfSpace Hddm_u~{i}~{iDom}~{iEdge}; }
+    For iSide In {0:1}
+      { Name du~{i}~{iDom}~{iEdge}~{iSide}; Type Local; NameOfSpace Hddm_du~{i}~{iDom}~{iEdge}~{iSide}; }
+    EndFor
+    EndFor
+    EndFor
+    For iCorner In {0:3}
+    iEdgeSide~{0} = (iCorner+3)%4;
+    iEdgeSide~{1} = iCorner;
+    For i In {1:nPade~{iDom}~{iEdgeSide~{0}}}
+    For j In {1:nPade~{iDom}~{iEdgeSide~{1}}}
+      { Name u~{i}~{j}~{iDom}~{iCorner}; Type Local; NameOfSpace Hddm_u~{i}~{j}~{iDom}~{iCorner}; }
+    EndFor
+    EndFor
+    EndFor
+    }
+    Equation {
+
+// Field on SURFACE
+
+      Galerkin { [ Dof{d u~{iDom}}                                , {d u~{iDom}} ]; In Omega~{iDom};         Jacobian JSur; Integration I1; }
+      Galerkin { [ - k[]^2 * Dof{u~{iDom}}                          , {u~{iDom}} ]; In Omega~{iDom};         Jacobian JSur; Integration I1; }
+      Galerkin { [ ($PhysicalSource ? -f_sou~{iDom} : 0)            , {u~{iDom}} ]; In GammaPoint~{iDom};    Jacobian JSur; Integration I1; }
+      Galerkin { [ ($PhysicalSource ? -f_neu[] : 0)                 , {u~{iDom}} ]; In GammaN~{iDom};        Jacobian JLin; Integration I1; }
+    For iEdge In {0:3}
+      Galerkin { [ - Dof{du~{iDom}~{iEdge}}                         , {u~{iDom}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JLin; Integration I1; }
+    If (isEdgeTBC~{iDom}~{iEdge})
+      Galerkin { [ ($ArtificialSource ? -g_in~{iDom}~{iEdge}[] : 0) , {u~{iDom}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JLin; Integration I1; }
+    EndIf
+    EndFor
+
+// DtN on EDGES
+
+    For iEdge In {0:3}
+      Galerkin { [   Dof{du~{iDom}~{iEdge}}                                                                                              , {du~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JLin; Integration I1; }
+      Galerkin { [ - paramCurv1~{iDom}~{iEdge}[] * Dof{u~{iDom}}                                                                         , {du~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JLin; Integration I1; }
+      Galerkin { [ - paramCurv2~{iDom}~{iEdge}[] * Dof{d u~{iDom}}                                                                     , {d du~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JLin; Integration I1; }
+      Galerkin { [ - I[]*k[]*alphaPade~{iDom}~{iEdge}[] * Dof{u~{iDom}}                                                                  , {du~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JLin; Integration I1; }
+    For i In {1:nPade~{iDom}~{iEdge}}
+      Galerkin { [ - I[]*k[]*alphaPade~{iDom}~{iEdge}[] * 2./mPade~{iDom}~{iEdge} * cPade~{i}~{iDom}~{iEdge} * Dof{u~{iDom}}             , {du~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JLin; Integration I1; }
+      Galerkin { [ - I[]*k[]*alphaPade~{iDom}~{iEdge}[] * 2./mPade~{iDom}~{iEdge} * cPade~{i}~{iDom}~{iEdge} * Dof{u~{i}~{iDom}~{iEdge}} , {du~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JLin; Integration I1; }
+    EndFor
+    EndFor
+
+// Fields on EDGES
+
+    For iEdge In {0:3}
+    iEdgeSide~{0} = (iEdge+3)%4;
+    iEdgeSide~{1} = (iEdge+1)%4;
+    iCornerSide~{0} = iEdge;
+    iCornerSide~{1} = (iEdge+1)%4;
+    iSideSide~{0} = 1;
+    iSideSide~{1} = 0;
+    For i In {1:nPade~{iDom}~{iEdge}}
+
+      Galerkin { [ Dof{d u~{i}~{iDom}~{iEdge}}                                                                 , {d u~{i}~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge};                Jacobian JLin; Integration I1; }
+      Galerkin { [ - kEps~{iDom}~{iEdge}[]^2*(alphaPade~{iDom}~{iEdge}[]^2*cPade~{i}~{iDom}~{iEdge}+1) * Dof{u~{i}~{iDom}~{iEdge}} , {u~{i}~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge};                Jacobian JLin; Integration I1; }
+      Galerkin { [ - kEps~{iDom}~{iEdge}[]^2*alphaPade~{iDom}~{iEdge}[]^2*(cPade~{i}~{iDom}~{iEdge}+1) * Dof{u~{iDom}}             , {u~{i}~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge};                Jacobian JLin; Integration I1; }
+    For iSide In {0:1}
+    If ((CROSSPOINT_TREAT == CROSSPOINT_Full) || (isEdgeABC~{iDom}~{iEdge} && isEdgeABC~{iDom}~{iEdgeSide~{iSide}}))
+      Galerkin { [ - Dof{du~{i}~{iDom}~{iEdge}~{iSide}}                                                          , {u~{i}~{iDom}~{iEdge}} ]; In Corner~{iDom}~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    If (isEdgeTBC~{iDom}~{iEdgeSide~{iSide}})
+      Galerkin { [ ($ArtificialSource ? -g_in~{i}~{iDom}~{iEdgeSide~{iSide}}~{iSideSide~{iSide}}[] : 0)          , {u~{i}~{iDom}~{iEdge}} ]; In Corner~{iDom}~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    EndIf
+    EndIf
+    EndFor
+    EndFor
+    EndFor
+
+// DtN at CORNERS
+
+    For iEdge In {0:3}
+    iEdgeSide~{0} = (iEdge+3)%4;
+    iEdgeSide~{1} = (iEdge+1)%4;
+    iCornerSide~{0} = iEdge;
+    iCornerSide~{1} = (iEdge+1)%4;
+    For i In {1:nPade~{iDom}~{iEdge}}
+
+    For iSide In {0:1}
+      Galerkin { [   Dof{du~{i}~{iDom}~{iEdge}~{iSide}}                                                                                                                                        , {du~{i}~{iDom}~{iEdge}~{iSide}} ]; In Corner~{iDom}~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+      Galerkin { [ - I[]*k[]*alphaPade~{iDom}~{iEdgeSide~{iSide}}[] * Dof{u~{i}~{iDom}~{iEdge}}                                                                                                , {du~{i}~{iDom}~{iEdge}~{iSide}} ]; In Corner~{iDom}~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    For j In {1:nPade~{iDom}~{iEdgeSide~{iSide}}}
+      Galerkin { [ - I[]*k[]*alphaPade~{iDom}~{iEdgeSide~{iSide}}[] * 2./mPade~{iDom}~{iEdgeSide~{iSide}} * cPade~{j}~{iDom}~{iEdgeSide~{iSide}} * Dof{u~{i}~{iDom}~{iEdge}}                   , {du~{i}~{iDom}~{iEdge}~{iSide}} ]; In Corner~{iDom}~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    If (iSide == 0)
+      Galerkin { [ - I[]*k[]*alphaPade~{iDom}~{iEdgeSide~{iSide}}[] * 2./mPade~{iDom}~{iEdgeSide~{iSide}} * cPade~{j}~{iDom}~{iEdgeSide~{iSide}} * Dof{u~{j}~{i}~{iDom}~{iCornerSide~{iSide}}} , {du~{i}~{iDom}~{iEdge}~{iSide}} ]; In Corner~{iDom}~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    Else
+      Galerkin { [ - I[]*k[]*alphaPade~{iDom}~{iEdgeSide~{iSide}}[] * 2./mPade~{iDom}~{iEdgeSide~{iSide}} * cPade~{j}~{iDom}~{iEdgeSide~{iSide}} * Dof{u~{i}~{j}~{iDom}~{iCornerSide~{iSide}}} , {du~{i}~{iDom}~{iEdge}~{iSide}} ]; In Corner~{iDom}~{iCornerSide~{iSide}}; Jacobian JVol; Integration I1; }
+    EndIf
+    EndFor
+    EndFor
+
+    EndFor
+    EndFor
+
+// Fields at CORNERS
+
+    For iCorner In {0:3}
+    iEdge0 = (iCorner+3)%4;
+    iEdge1 = iCorner;
+    For i In {1:nPade~{iDom}~{iEdge0}}
+    For j In {1:nPade~{iDom}~{iEdge1}}
+      Galerkin { [ Dof{u~{i}~{j}~{iDom}~{iCorner}}                                                                                                                                                                                , {u~{i}~{j}~{iDom}~{iCorner}} ]; In Corner~{iDom}~{iCorner}; Jacobian JVol; Integration I1; }
+      Galerkin { [ alphaPade~{iDom}~{iEdge1}[]^2*(cPade~{j}~{iDom}~{iEdge1}+1) * Dof{u~{i}~{iDom}~{iEdge0}} / (alphaPade~{iDom}~{iEdge0}[]^2*cPade~{i}~{iDom}~{iEdge0} + alphaPade~{iDom}~{iEdge1}[]^2*cPade~{j}~{iDom}~{iEdge1} + 1) , {u~{i}~{j}~{iDom}~{iCorner}} ]; In Corner~{iDom}~{iCorner}; Jacobian JVol; Integration I1; }
+      Galerkin { [ alphaPade~{iDom}~{iEdge0}[]^2*(cPade~{i}~{iDom}~{iEdge0}+1) * Dof{u~{j}~{iDom}~{iEdge1}} / (alphaPade~{iDom}~{iEdge0}[]^2*cPade~{i}~{iDom}~{iEdge0} + alphaPade~{iDom}~{iEdge1}[]^2*cPade~{j}~{iDom}~{iEdge1} + 1) , {u~{i}~{j}~{iDom}~{iCorner}} ]; In Corner~{iDom}~{iCorner}; Jacobian JVol; Integration I1; }
+    EndFor
+    EndFor
+    EndFor
+
+    }
+  }
+
+EndFor // loop on iDom
+
+// ====================================================================================================
+// SOLVER for outgoing data for side 'iEdge' of 'iDom'
+// ====================================================================================================
+
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+For i1 In {0:#ListOfEdgesWithTBC~{iDom}()-1}
+iEdge = ListOfEdgesWithTBC~{iDom}(i1);
+iEdgeSide~{0} = (iEdge+3)%4;
+iEdgeSide~{1} = (iEdge+1)%4;
+iCornerSide~{0} = iEdge;
+iCornerSide~{1} = (iEdge+1)%4;
+iSideSide~{0} = 1;
+iSideSide~{1} = 0;
+
+  { Name FormSur~{iDom}~{iEdge}; Type FemEquation;
+    Quantity {
+      { Name du~{iDom}~{iEdge};    Type Local; NameOfSpace Hddm_du~{iDom}~{iEdge}; }
+      { Name g_out~{iDom}~{iEdge}; Type Local; NameOfSpace Hddm_g_out~{iDom}~{iEdge}; }
+    For iSide In {0:1}
+    For i In {1:nPade~{iDom}~{iEdgeSide~{iSide}}}
+      { Name du~{i}~{iDom}~{iEdgeSide~{iSide}}~{iSideSide~{iSide}}; Type Local; NameOfSpace Hddm_du~{i}~{iDom}~{iEdgeSide~{iSide}}~{iSideSide~{iSide}}; }
+      { Name g_out~{i}~{iDom}~{iEdge}~{iSide};                      Type Local; NameOfSpace Hddm_g_out~{i}~{iDom}~{iEdge}~{iSide}; }
+    EndFor
+    EndFor
+    }
+    Equation {
+
+      Galerkin { [ Dof{g_out~{iDom}~{iEdge}}                     , {g_out~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JSur; Integration I1; }
+      Galerkin { [ $ArtificialSource ? g_in~{iDom}~{iEdge}[] : 0 , {g_out~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JSur; Integration I1; }
+      Galerkin { [ 2 * {du~{iDom}~{iEdge}}                       , {g_out~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JSur; Integration I1; }
+
+    For iSide In {0:1}
+    If (CROSSPOINT_TREAT == CROSSPOINT_Full)
+    For i In {1:nPade~{iDom}~{iEdgeSide~{iSide}}}
+      Galerkin { [ Dof{g_out~{i}~{iDom}~{iEdge}~{iSide}}                       , {g_out~{i}~{iDom}~{iEdge}~{iSide}} ]; In Corner~{iDom}~{iCornerSide~{iSide}}; Jacobian JLin; Integration I1; }
+      Galerkin { [ $ArtificialSource ? g_in~{i}~{iDom}~{iEdge}~{iSide}[] : 0   , {g_out~{i}~{iDom}~{iEdge}~{iSide}} ]; In Corner~{iDom}~{iCornerSide~{iSide}}; Jacobian JLin; Integration I1; }
+      Galerkin { [ 2 * {du~{i}~{iDom}~{iEdgeSide~{iSide}}~{iSideSide~{iSide}}} , {g_out~{i}~{iDom}~{iEdge}~{iSide}} ]; In Corner~{iDom}~{iCornerSide~{iSide}}; Jacobian JLin; Integration I1; }
+    EndFor
+    EndIf
+    EndFor
+
+    }
+  }
+
+EndFor // loop on iEdge
+EndFor // loop on iDom
+
+}
+
+// ====================================================================================================
+// RESOLUTION
+// ====================================================================================================
+
+Resolution {
+  { Name ResoRef;
+    System {
+        { Name SysRef; NameOfFormulation FormRef; Type Complex; NameOfMesh Sprintf[StrCat[MSH_NAME, ".msh"]]; }
+    }
+    Operation {
+        Generate[SysRef]; Solve[SysRef]; PostOperation[PostOpRef];
+    }
+  }
+
+  { Name ResoDDM;
+    System {
+      If(FLAG_ERROR == ERROR_NUM)
+        { Name SysRef; NameOfFormulation FormRef; Type Complex; NameOfMesh Sprintf[StrCat[MSH_NAME, ".msh"]]; }
+      EndIf
+      For i0 In {0:#ListOfSubdom()-1}
+      iDom = ListOfSubdom(i0);
+        { Name SysVol~{iDom}; NameOfFormulation FormVol~{iDom}; Type Complex; NameOfMesh Sprintf[StrCat[MSH_NAME, "_%g.msh"], iDom]; }
+      For i1 In {0:#ListOfEdgesWithTBC~{iDom}()-1}
+      iEdge = ListOfEdgesWithTBC~{iDom}(i1);
+        { Name SysSur~{iDom}~{iEdge}; NameOfFormulation FormSur~{iDom}~{iEdge}; Type Complex; NameOfMesh Sprintf[StrCat[MSH_NAME, "_%g.msh"], iDom]; }
+      EndFor
+      EndFor
+    }
+    Operation {
+
+      If(FLAG_ERROR == ERROR_NUM)
+        Generate[SysRef]; Solve[SysRef]; PostOperation[PostOpRef];
+      EndIf
+
+//    For MAXITeff In {1:50}
+
+      // =========================================================
+      // BUILD local part of distributed rhs b for Krylov solver using physical sources only, and update surface data
+      // =========================================================
+
+      SetCommSelf;
+
+      // SolveVolumePDE (without artificial/with physical sources)
+      Evaluate[$ArtificialSource = 0];
+      Evaluate[$PhysicalSource = 1];
+      For i0 In {0:#ListOfSubdom()-1}
+      iDom = ListOfSubdom(i0);
+        UpdateConstraint[SysVol~{iDom}, GammaD~{iDom}, Assign];
+        Generate[SysVol~{iDom}];
+        Solve[SysVol~{iDom}];
+      EndFor
+
+    If(#ListOfSubdom() > 1)
+
+      // SolveSurfacePDE
+      For i0 In {0:#ListOfSubdom()-1}
+      iDom = ListOfSubdom(i0);
+      For i1 In {0:#ListOfEdgesWithTBC~{iDom}()-1}
+      iEdge = ListOfEdgesWithTBC~{iDom}(i1);
+        Generate[SysSur~{iDom}~{iEdge}];
+        Solve[SysSur~{iDom}~{iEdge}];
+      EndFor
+      EndFor
+
+      // UpdateSurfaceFields
+      For i0 In {0:#ListOfSubdom()-1}
+      iDom = ListOfSubdom(i0);
+      For i1 In {0:#ListOfEdgesWithTBC~{iDom}()-1}
+      iEdge = ListOfEdgesWithTBC~{iDom}(i1);
+        PostOperation[PostOpSur~{iDom}~{iEdge}];
+      EndFor
+      EndFor
+
+      // =========================================================
+      // Krylov solver using artificial sources only. IterativeLinearSolver solves (I-A) g = b
+      // ListOfFields() initially stores the local part of b; then stores each local part of iterate g^n.
+      // =========================================================
+
+      // UpdateConstraint (with artificial/without physical sources)
+      Evaluate[$ArtificialSource = 1];
+      Evaluate[$PhysicalSource = 0];
+      For i0 In {0:#ListOfSubdom()-1}
+      iDom = ListOfSubdom(i0);
+        UpdateConstraint[SysVol~{iDom}, GammaD~{iDom}, Assign];
+      EndFor
+
+      SetCommWorld;
+
+      IterativeLinearSolver["I-A", SOLVER, TOL, MAXIT, RESTART, {ListOfFields()}, {ListOfConnectedFields()}, {}]
+//      IterativeLinearSolver["I-A", SOLVER, TOL, MAXITeff, RESTART, {ListOfFields()}, {ListOfConnectedFields()}, {}]
+      {
+
+        SetCommSelf;
+
+        // SolveVolumePDE
+        For i0 In {0:#ListOfSubdom()-1}
+        iDom = ListOfSubdom(i0);
+          GenerateRHSGroup[SysVol~{iDom}, Region[{GammaD~{iDom}, Sigma~{iDom}, Corner~{iDom}}] ];
+          SolveAgain[SysVol~{iDom}];
+        EndFor
+
+        // SolveSurfacePDE
+        For i0 In {0:#ListOfSubdom()-1}
+        iDom = ListOfSubdom(i0);
+        For i1 In {0:#ListOfEdgesWithTBC~{iDom}()-1}
+        iEdge = ListOfEdgesWithTBC~{iDom}(i1);
+          GenerateRHSGroup[SysSur~{iDom}~{iEdge}, Region[{Sigma~{iDom}~{iEdge}, SigmaSide~{iDom}~{iEdge}~{0}, SigmaSide~{iDom}~{iEdge}~{1}}]];
+          SolveAgain[SysSur~{iDom}~{iEdge}];
+        EndFor
+        EndFor
+
+        // UpdateSurfaceFields
+        For i0 In {0:#ListOfSubdom()-1}
+        iDom = ListOfSubdom(i0);
+        For i1 In {0:#ListOfEdgesWithTBC~{iDom}()-1}
+        iEdge = ListOfEdgesWithTBC~{iDom}(i1);
+          PostOperation[PostOpSur~{iDom}~{iEdge}];
+        EndFor
+        EndFor
+
+        SetCommWorld;
+
+        Barrier;
+        Test[ $KSPIteration == 0 ]{
+          Print[ {$KSPResidual}, Format "%12.5e", File Sprintf("out/STAT_residualInit_%g_%g_%g_%g_%g_%g.dat", FLAG_PBM, N_DOM, nPadeABC, thetaPadeInputABC, nPadeTBC, CROSSPOINT_TREAT) ];
+          //Print[ {$KSPResidual}, Format "%12.5e", File Sprintf("out/STAT_residualInit_%g_%g_%g_%g_%g_%g_%g.dat", FLAG_PBM, N_DOM, nPadeABC, thetaPadeInputABC, nPadeTBC, thetaPadeInputTBC, CROSSPOINT_TREAT) ];
+        }
+        //Print[ {$KSPResidual}, Format "%12.5e", File Sprintf("out/STAT_solverRes_%g_%g_%g_%g_%g_%g_%g_%g_%g_%g.dat", FLAG_PBM, N_DOM, nPadeABC, thetaPadeInputABC, nPadeTBC, thetaPadeInputTBC, CROSSPOINT_TREAT, WAVENUMBER, N_LAMBDAtextsf) ];
+      }
+
+      //Print[ {FLAG_PBM, N_DOM, nPadeABC, thetaPadeInputABC, nPadeTBC, thetaPadeInputTBC, CROSSPOINT_TREAT, WAVENUMBER, N_LAMBDA, $KSPIts}, Format "%g %g %g %g %g %g %g %g %g %g", File Sprintf("out/STAT_solverNumIter.dat") ];
+      Evaluate[$NumIter = $KSPIts];
+      Evaluate[$ResidualFinal = $KSPResidual];
+
+      //Print[ {$KSPResidual}, Format "%12.5e", File Sprintf("out/STAT_solverRes_%g_%g_%g_%g_%g_%g_%g_%g_%g.dat", FLAG_PBM, N_DOM, nPadeABC, thetaPadeInputABC, nPadeTBC, thetaPadeInputTBC, CROSSPOINT_TREAT, WAVENUMBER, N_LAMBDA) ];
+
+    EndIf
+
+      SetCommSelf;
+
+      // Compute final solution
+      Evaluate[$ArtificialSource = 1];
+      Evaluate[$PhysicalSource = 1];
+      For i0 In {0:#ListOfSubdom()-1}
+      iDom = ListOfSubdom(i0);
+        UpdateConstraint[SysVol~{iDom}, GammaD~{iDom}, Assign];
+        Generate[SysVol~{iDom}];
+        Solve[SysVol~{iDom}];
+      EndFor
+
+      // Save final solution
+      For i0 In {0:#ListOfSubdom()-1}
+      iDom = ListOfSubdom(i0);
+        PostOperation[PostOpVol~{iDom}];
+      EndFor
+
+      SetCommWorld;
+
+      // Compute error
+      If(FLAG_ERROR != ERROR_NO)
+        Evaluate[$error2Var = 0];
+        Evaluate[$energy2Var = 0];
+        For i0 In {0:#ListOfSubdom()-1}
+        iDom = ListOfSubdom(i0);
+          PostOperation[PostOpStat~{iDom}];
+          Evaluate[$error2Var = $error2Var + $error2Var~{iDom}];
+          Evaluate[$energy2Var = $energy2Var + $energy2Var~{iDom}];
+        EndFor
+      EndIf
+      PostOperation[PostOpStat];
+
+//    EndFor
+
+    }
+  }
+}
+
+
+// ====================================================================================================
+// POST-PROCESSING
+// ====================================================================================================
+
+PostProcessing {
+
+// REFERENCE SOLVER
+
+  { Name PostProRef; NameOfFormulation FormRef;
+    PostQuantity {
+//      { Name cNumREF; Value { Local { [ cData[] ];              In Omega; Jacobian JVol; }}}
+      { Name uNumREF; Value { Local { [ {uRef} ];               In Omega; Jacobian JVol; }}}
+    If((FLAG_PBM == PBM_SCATT_CIRC) || (FLAG_PBM == PBM_SCATT_RECT))
+      { Name uRefREF; Value { Local { [ f_ref[] ];              In Omega; Jacobian JVol; }}}
+      { Name uErrREF; Value { Local { [ Norm[f_ref[]-{uRef}] ]; In Omega; Jacobian JVol; }}}
+    EndIf
+    }
+  }
+
+// SUBDOMAIN SOLVER
+
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+  { Name PostProVol~{iDom}; NameOfFormulation FormVol~{iDom};
+    PostQuantity {
+      { Name uNum~{iDom};    Value { Local { [ {u~{iDom}} ];                    In Omega~{iDom}; Jacobian JVol; }}}
+    If(FLAG_ERROR != ERROR_NO)
+      { Name uRef~{iDom};    Value { Local { [ f_refErr[] ];                    In Omega~{iDom}; Jacobian JVol; }}}
+      { Name uErr~{iDom};    Value { Local { [ Norm[f_refErr[]-{u~{iDom}}]^2 ]; In Omega~{iDom}; Jacobian JVol; }}}
+    EndIf
+    }
+  }
+EndFor
+
+// INTERFACE CONDITION
+
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+For i1 In {0:#ListOfEdgesWithTBC~{iDom}()-1}
+iEdge = ListOfEdgesWithTBC~{iDom}(i1);
+iEdgeSide~{0} = (iEdge+3)%4;
+iEdgeSide~{1} = (iEdge+1)%4;
+  { Name PostProSur~{iDom}~{iEdge}; NameOfFormulation FormSur~{iDom}~{iEdge};
+    PostQuantity {
+      { Name g_out~{iDom}~{iEdge}; Value { Local { [ {g_out~{iDom}~{iEdge}} ]; In Sigma~{iDom}~{iEdge}; Jacobian JSur; }}}
+    For iSide In {0:1}
+    If (CROSSPOINT_TREAT == CROSSPOINT_Full)
+    For i In {1:nPade~{iDom}~{iEdgeSide~{iSide}}}
+      { Name g_out~{i}~{iDom}~{iEdge}~{iSide}; Value { Local { [ {g_out~{i}~{iDom}~{iEdge}~{iSide}} ]; In SigmaSide~{iDom}~{iEdge}~{iSide}; Jacobian JLin; }}}
+    EndFor
+    EndIf
+    EndFor
+    }
+  }
+EndFor
+EndFor
+
+// STATS (subdomain)
+
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+  { Name PostProStat~{iDom}; NameOfFormulation FormVol~{iDom};
+    PostQuantity {
+    If(FLAG_ERROR != ERROR_NO)
+      { Name energy2~{iDom}; Value { Integral { [ Norm[f_refErr[]]^2 ];            In Omega~{iDom}; Jacobian JVol; Integration I1; }}}
+      { Name error2~{iDom};  Value { Integral { [ Norm[f_refErr[]-{u~{iDom}}]^2 ]; In Omega~{iDom}; Jacobian JVol; Integration I1; }}}
+    EndIf
+    }
+  }
+EndFor
+
+// STATS (global)
+
+If (MPI_Rank == 0)
+  { Name PostProStat; NameOfFormulation FormVol~{0};
+    PostQuantity {
+    If(FLAG_ERROR != ERROR_NO)
+      { Name energAbs; Value { Term { Type Global; [Sqrt[$energy2Var]] ;            In Omega~{0}; Jacobian JVol; }}}
+      { Name errorAbs; Value { Term { Type Global; [Sqrt[$error2Var]] ;             In Omega~{0}; Jacobian JVol; }}}
+      { Name errorRel; Value { Term { Type Global; [Sqrt[$error2Var/$energy2Var]] ; In Omega~{0}; Jacobian JVol; }}}
+    EndIf
+      { Name numIter;  Value { Term { Type Global; [$NumIter];                      In Omega~{0}; Jacobian JVol; }}}
+      { Name residual; Value { Term { Type Global; [$ResidualFinal];                In Omega~{0}; Jacobian JVol; }}}
+    }
+  }
+EndIf
+}
+
+// ====================================================================================================
+// POST-OPERATIONS
+// ====================================================================================================
+
+PostOperation {
+
+// REFERENCE SOLVER
+
+  { Name PostOpRef; NameOfPostProcessing PostProRef;
+    Operation {
+//      Print [uNumREF, OnElementsOf Omega, StoreInField (FIELD_REF), AtGaussPoints NumGaussPointsTRI ];
+      Print [uNumREF, OnElementsOf Omega, StoreInField (FIELD_REF)];
+//      Print [uNumREF, OnElementsOf Omega, File StrCat(DIR, "uNumREF.pos") ];
+//      Print [cNumREF, OnElementsOf Omega, File StrCat(DIR, "cNumREF.pos") ];
+//    If((FLAG_PBM == PBM_SCATT_CIRC) || (FLAG_PBM == PBM_SCATT_RECT))
+//      Print [uRefREF, OnElementsOf Omega, File StrCat(DIR, "uRef.pos") ];
+//      Print [uErrREF, OnElementsOf Omega, File StrCat(DIR, "uErr.pos") ];
+//    EndIf
+    }
+  }
+
+// SUB-DOMAIN SOLVER
+
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+  { Name PostOpVol~{iDom}; NameOfPostProcessing PostProVol~{iDom};
+    Operation {
+      Print [uNum~{iDom}, OnElementsOf Omega~{iDom}, File StrCat(DIR, Sprintf("uNum_%g.pos", iDom))];
+    If(FLAG_ERROR != ERROR_NO)
+      Print [uRef~{iDom}, OnElementsOf Omega~{iDom}, File StrCat(DIR, Sprintf("uRef_%g.pos", iDom))];
+      Print [uErr~{iDom}, OnElementsOf Omega~{iDom}, File StrCat(DIR, Sprintf("uErr_%g.pos", iDom))];
+    EndIf
+    }
+  }
+EndFor
+
+// INTERFACE CONDITION
+
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+For i1 In {0:#ListOfEdgesWithTBC~{iDom}()-1}
+iEdge = ListOfEdgesWithTBC~{iDom}(i1);
+iEdgeSide~{0} = (iEdge+3)%4;
+iEdgeSide~{1} = (iEdge+1)%4;
+  { Name PostOpSur~{iDom}~{iEdge}; NameOfPostProcessing PostProSur~{iDom}~{iEdge};
+    Operation {
+      Print [ g_out~{iDom}~{iEdge}, OnElementsOf Sigma~{iDom}~{iEdge}, StoreInField (FIELD~{iDom}~{iEdge}), AtGaussPoints NumGaussPointsLIN ];
+    For iSide In {0:1}
+    If (CROSSPOINT_TREAT == CROSSPOINT_Full)
+    For i In {1:nPade~{iDom}~{iEdgeSide~{iSide}}}
+      Print [ g_out~{i}~{iDom}~{iEdge}~{iSide}, OnElementsOf SigmaSide~{iDom}~{iEdge}~{iSide}, StoreInField (FIELD~{i}~{iDom}~{iEdge}~{iSide}), AtGaussPoints NumGaussPointsPNT ];
+    EndFor
+    EndIf
+    EndFor
+    }
+  }
+EndFor
+EndFor
+
+// STATS (subdomain)
+
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+  { Name PostOpStat~{iDom}; NameOfPostProcessing PostProStat~{iDom};
+    Operation {
+    If(FLAG_ERROR != ERROR_NO)
+      Print [energy2~{iDom}[Omega~{iDom}], OnRegion Omega~{iDom}, Format Table, StoreInVariable $energy2Var~{iDom}];
+      Print [error2~{iDom}[Omega~{iDom}],  OnRegion Omega~{iDom}, Format Table, StoreInVariable $error2Var~{iDom}];
+    EndIf
+    }
+  }
+EndFor
+
+// STATS (global)
+
+If (MPI_Rank == 0)
+  { Name PostOpStat; NameOfPostProcessing PostProStat;
+    Operation {
+      tmp0 = Sprintf("out/STAT_energAbs_%g_%g_%g_%g_%g_%g.dat",      FLAG_PBM, N_DOM, nPadeABC, thetaPadeInputABC, nPadeTBC, CROSSPOINT_TREAT);
+      tmp1 = Sprintf("out/STAT_errorAbs_%g_%g_%g_%g_%g_%g.dat",      FLAG_PBM, N_DOM, nPadeABC, thetaPadeInputABC, nPadeTBC, CROSSPOINT_TREAT);
+      tmp2 = Sprintf("out/STAT_errorRel_%g_%g_%g_%g_%g_%g.dat",      FLAG_PBM, N_DOM, nPadeABC, thetaPadeInputABC, nPadeTBC, CROSSPOINT_TREAT);
+      tmp3 = Sprintf("out/STAT_numIter_%g_%g_%g_%g_%g_%g.dat",       FLAG_PBM, N_DOM, nPadeABC, thetaPadeInputABC, nPadeTBC, CROSSPOINT_TREAT);
+      tmp4 = Sprintf("out/STAT_residualFinal_%g_%g_%g_%g_%g_%g.dat", FLAG_PBM, N_DOM, nPadeABC, thetaPadeInputABC, nPadeTBC, CROSSPOINT_TREAT);
+    If(FLAG_ERROR != ERROR_NO)
+      Print [energAbs, OnRegion Omega~{0}, SendToServer "Output/1energAbs",       Format Table, File > tmp0];
+      Print [errorAbs, OnRegion Omega~{0}, SendToServer "Output/2errorAbs",       Format Table, File > tmp1];
+      Print [errorRel, OnRegion Omega~{0}, SendToServer "Output/3errorRel",       Format Table, File > tmp2];
+    EndIf
+      Print [numIter,  OnRegion Omega~{0}, SendToServer "Output/4Num iterations", Format Table, File > tmp3];
+      Print [residual, OnRegion Omega~{0}, SendToServer "Output/5Residual",       Format Table, File > tmp4];
+    }
+  }
+EndIf
+}
+
+// ====================================================================================================
+
+DefineConstant[
+  R_ = {"ResoDDM", Choices {"ResoRef", "ResoDDM"}, Name "GetDP/1ResolutionChoices", Visible 1},
+  C_ = {"-solve -v 3 -bin -v2 -ksp_monitor", Name "GetDP/9ComputeCommand", Visible 0}
+  P_ = {"", Name "GetDP/2PostOperationChoices", Visible 1}
+];

HelmholtzDDMwithCrosspoints/main.dat

+DefineConstant[
+  MSH_BASE_NAME = "mesh",
+  DIR = "out/"
+];
+MSH_NAME = StrCat(DIR, MSH_BASE_NAME);
+
+// =================================================================================================
+// SPECIFIC PROBLEM
+// =================================================================================================
+
+PBM_SCATT_CIRC = 1;
+PBM_SCATT_RECT = 2;
+PBM_MARMOUSI   = 3;
+
+DefineConstant[
+  FLAG_PBM = {PBM_SCATT_RECT, Name "Input/01Benchmark", Highlight "Blue",
+              GmshOption "Reset", Autocheck 0,
+              Choices {PBM_SCATT_CIRC = "Scattering in circular domain",
+                       PBM_SCATT_RECT = "Scattering in rectangular domain",
+                       PBM_MARMOUSI   = "Benchmark Marmousi"} }
+];
+
+If (FLAG_PBM == PBM_SCATT_CIRC)
+  LinkGeo = "circularDom.geo";
+  LinkPro = "circularDom.pro";
+EndIf
+If ((FLAG_PBM == PBM_SCATT_RECT) ||  (FLAG_PBM == PBM_MARMOUSI))
+  LinkGeo = "rectangularDom.geo";
+  LinkPro = "rectangularDom.pro";
+EndIf
+LinkDDM = "ddmHelmholtzCrossPoints.pro";
+
+ERROR_NO  = 0;
+ERROR_ANA = 1;
+ERROR_NUM = 2;
+
+DefineConstant[
+  FLAG_ERROR = {ERROR_NUM, Name "Input/02Error", Autocheck 1,
+              Choices {ERROR_NO  = "No",
+                       ERROR_ANA = "... vs analytic reference",
+                       ERROR_NUM = "... vs numerical reference"} },
+  ORDER = {1,
+    Name "Input/05Polynomial order",
+    Choices {1 = "First-order", 2 = "Second-order"}}
+];

HelmholtzDDMwithCrosspoints/main.geo

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

HelmholtzDDMwithCrosspoints/main.pro

+Include "main.dat" ;
+
+// =================================================================================================
+// GENERAL OPTIONS
+// =================================================================================================
+
+DefineConstant[
+  SOLVER = {"gmres", Choices {"gmres", "fgmres", "bcgs", "print", "jacobi"},
+    Name "Iterative Solver/0Solver"},
+  TOLlog10 = {-6, Max -1, Min -16, Step -1,
+    Name "Iterative Solver/1Tolerance (log10)"},
+  TOL = 10^(TOLlog10),
+  MAXIT = {500, Min 1, Step 1, Max 100000,
+    Name "Iterative Solver/2Max. iterations"},
+  RESTART_MAXIT = {0, Choices {0,1}, Visible 0,
+    Name "Iterative Solver/31Force Restart = Max. iterations"},
+  RESTART = {RESTART_MAXIT ? MAXIT : MAXIT, Min 0, Max 100000, Step 1, Visible 0,
+    Name "Iterative Solver/30Restart", ReadOnly RESTART_MAXIT }
+];
+
+If(RESTART > MAXIT || RESTART == 0)
+  RESTART = MAXIT;
+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; }}}
+}
+
+If (ORDER==1)
+  NumGaussPointsPNT = 1;
+  NumGaussPointsLIN = 4; // NumberOfPoints
+  NumGaussPointsTRI = 6; // NumberOfPoints
+EndIf
+If (ORDER==2)
+  NumGaussPointsPNT = 1;
+  NumGaussPointsLIN = 6;  // NumberOfPoints
+  NumGaussPointsTRI = 12; // NumberOfPoints
+EndIf
+
+Integration {
+  { Name I1;
+    Case {
+      { Type Gauss;
+        Case {
+          { GeoElement Point;      NumberOfPoints NumGaussPointsPNT; }
+          { GeoElement Line;       NumberOfPoints NumGaussPointsLIN; }
+          { GeoElement Line2;      NumberOfPoints NumGaussPointsLIN; }
+          { GeoElement Triangle;   NumberOfPoints NumGaussPointsTRI; }
+          { GeoElement Triangle2;  NumberOfPoints NumGaussPointsTRI; } // 1 3 4 6 7 12 13 16
+          { GeoElement Quadrangle; NumberOfPoints 7; }
+        }
+      }
+    }
+  }
+}
+
+// =================================================================================================
+// SPECIFIC
+// =================================================================================================
+
+Include Str[LinkPro];
+Include Str[LinkDDM];

HelmholtzDDMwithCrosspoints/rectangularDom.dat

+// =================================================================================================
+
+// Geometry: Scattering case
+If(FLAG_PBM == PBM_SCATT_RECT)
+
+// Subdomains
+DefineConstant[
+  Lx = {7.5, Min 0.01, Max 20, Step 0.01, Name "Input/1Model/40Domain length (Lx)"},
+  Ly = {7.5, Min 0.01, Max 20, Step 0.01, Name "Input/1Model/41Domain length (Ly)"},
+  X_SCA = { 1.25, Min -10,  Max 10, Step 0.01, Name "Input/1Model/42Scatterer position (x)"},
+  Y_SCA = { 1.25, Min -10,  Max 10, Step 0.01, Name "Input/1Model/43Scatterer position (y)"},
+  R_SCA = { 1.,   Min 0.01, Max 20, Step 0.01, Name "Input/1Model/44Scatterer radius"},
+  Nx_DOM  = {3, Min 1, Max 10, Step 1, Name "Input/1Model/50Nx subdomains"},
+  Ny_DOM  = {3, Min 1, Max 10, Step 1, Name "Input/1Model/51Ny subdomains"},
+  TWIST   = {0, Min 0, Max 1.5, Step 0.5, Name "Input/1Model/52 TWIST"}
+];
+
+  X0_DOM = -X_SCA;
+  Y0_DOM = -Y_SCA;
+//  X0_DOM = -Lx/Nx_DOM/2.;
+//  Y0_DOM = -Ly/Ny_DOM/2.;
+
+DefineConstant[
+//  LAMBDA     = {0.5, Min 0.1, Max 30, Step 0.1, Name "Input/03Wavelength"},
+  WAVENUMBER = {4*Pi, Min 1, Max 24, Step 1, Name "Input/03Wavenumber"},
+  N_LAMBDA   = {10,   Min 5, Max 30, Step 1, Name "Input/04Points per wavelength"}
+];
+
+  //WAVENUMBER = (2*Pi)/LAMBDA;
+  LAMBDA = (2*Pi)/WAVENUMBER;
+  FREQ = WAVENUMBER/(2*Pi);
+  LC = LAMBDA/N_LAMBDA;
+EndIf
+
+// =================================================================================================
+
+// Geometry: Marmousi case
+If(FLAG_PBM == PBM_MARMOUSI)
+
+// Subdomains
+DefineConstant[
+  Nx_DOM  = {15, Min 1, Max 10, Step 1, Name "Input/1Model/50Nx subdomains"},
+  Ny_DOM  = {4, Min 1, Max 10, Step 1, Name "Input/1Model/51Ny subdomains"}
+];
+  TWIST = 0;
+
+  Lx = 9192.;
+  Ly = 2904.;
+  X0_DOM = 0.;
+  Y0_DOM = -Ly;
+  R_SCA = 1000;
+
+  X_SOU = 4585.;
+  Y_SOU = -10.;
+  Lx_SDOM = Lx/Nx_DOM;
+  Ly_SDOM = Ly/Ny_DOM;
+  I_SOU = Ceil[(X_SOU-X0_DOM)/Lx_SDOM]-1;
+  J_SOU = Ceil[(Y_SOU-Y0_DOM)/Ly_SDOM]-1;
+  Idom_SOU = J_SOU*Nx_DOM + I_SOU;
+  //Printf("I_SOU: %g",I_SOU);
+  //Printf("J_SOU: %g",J_SOU);
+  //Printf("Idom_SOU: %g",Idom_SOU);
+
+DefineConstant[
+  FREQ     = {10, Min 0.1, Max 30, Step 0.1, Name "Input/03Frequency"}, //10
+  N_LAMBDA = {5,                             Name "Input/04Points per wavelength"} //10
+];
+
+  // meshing for shortest wavelength
+  cMin = 1500.;
+  WAVENUMBER = 2*Pi*FREQ/cMin;
+  LAMBDA = 2*Pi/WAVENUMBER ;
+  LC = LAMBDA/N_LAMBDA;
+EndIf
+
+// =================================================================================================
+
+N_DOM = Nx_DOM*Ny_DOM;
+
+// =================================================================================================
+
+If(FLAG_PBM == PBM_SCATT_RECT)
+  TAG_SCA = 1000;
+EndIf
+If(FLAG_PBM == PBM_MARMOUSI)
+  TAG_SOU = 1000;
+EndIf
+
+For iDom In {0:N_DOM-1}
+  TAG_DOM~{iDom} = 2000 + iDom;
+For iEdge In {0:3}
+  TAG_BND~{iDom}~{iEdge} = 3000 + 4*iDom + iEdge;
+EndFor
+For iCorner In {0:3}
+  TAG_CRN~{iDom}~{iCorner} = 4000 + 4*iDom + iCorner;
+EndFor
+EndFor

HelmholtzDDMwithCrosspoints/rectangularDom.geo

+Include "rectangularDom.dat";
+
+Mesh.CharacteristicLengthMin = LC;
+Mesh.CharacteristicLengthMax = LC;
+
+// =================================================================================================
+// Build SCATTERER or SOURCE
+// =================================================================================================
+
+If(FLAG_PBM == PBM_SCATT_RECT)
+  p0 = newp; Point(p0) = {     0,     0, 0};
+  p1 = newp; Point(p1) = { R_SCA,     0, 0};
+  p2 = newp; Point(p2) = {     0, R_SCA, 0};
+  p3 = newp; Point(p3) = {-R_SCA,     0, 0};
+  p4 = newp; Point(p4) = {     0,-R_SCA, 0};
+  c1 = newl; Circle(c1) = {p1, p0, p2};
+  c2 = newl; Circle(c2) = {p2, p0, p3};
+  c3 = newl; Circle(c3) = {p3, p0, p4};
+  c4 = newl; Circle(c4) = {p4, p0, p1};
+  llSca = newll; Line loop(llSca) = {-c1, -c2, -c3, -c4};
+EndIf
+
+If(FLAG_PBM == PBM_MARMOUSI)
+  p0 = newp; Point(p0) = {X_SOU, Y_SOU, 0};
+  Physical Point(TAG_SOU) = p0;
+EndIf
+
+// =================================================================================================
+// Build POINTS, STRAIGHT LINES, RECTANGLES
+// =================================================================================================
+
+For i In {0:Nx_DOM}
+For j In {0:Ny_DOM}
+  p~{i}~{j} = newp;
+  If((i == 1) && (j == 1))
+    xAdd = TWIST;
+    yAdd = 0;
+  ElseIf((i == 1) && (j == 2))
+    xAdd = 0;
+    yAdd = TWIST;
+  ElseIf((i == 2) && (j == 2))
+    xAdd = -TWIST;
+    yAdd = 0;
+  ElseIf((i == 2) && (j == 1))
+    xAdd = 0;
+    yAdd = -TWIST;
+  Else
+    xAdd = 0;
+    yAdd = 0;
+  EndIf
+  Point(p~{i}~{j}) = {X0_DOM + Lx*i/Nx_DOM + xAdd, Y0_DOM + Ly*j/Ny_DOM + yAdd, 0};
+EndFor
+EndFor
+
+For i In {0:Nx_DOM}
+For j In {0:Ny_DOM}
+  If (i < Nx_DOM)
+    lx~{i}~{j} = newl;
+    Line(lx~{i}~{j}) = {p~{i}~{j}, p~{i+1}~{j}};
+  EndIf
+  If (j < Ny_DOM)
+    ly~{i}~{j} = newl;
+    Line(ly~{i}~{j}) = {p~{i}~{j}, p~{i}~{j+1}};
+  EndIf
+EndFor
+EndFor
+
+For i In {0:Nx_DOM-1}
+For j In {0:Ny_DOM-1}
+  ll~{i}~{j} = newll;
+  Line loop(ll~{i}~{j}) = {lx~{i}~{j}, ly~{i+1}~{j}, -lx~{i}~{j+1}, -ly~{i}~{j}};
+  s~{i}~{j} = news;
+  If ((i == 0) && (j == 0) && (FLAG_PBM == PBM_SCATT_RECT))
+    Plane Surface(s~{i}~{j}) = {ll~{i}~{j}, llSca};
+  Else
+    Plane Surface(s~{i}~{j}) = {ll~{i}~{j}};
+  EndIf
+  If (FLAG_PBM == PBM_MARMOUSI)
+  If ((i == I_SOU) && (j == J_SOU))
+    Point {p0} In Surface {s~{i}~{j}};
+  EndIf
+  EndIf
+EndFor
+EndFor
+
+// =================================================================================================
+// Generate MESH for each subdomain and reference domain
+// =================================================================================================
+
+If(StrCmp(OnelabAction, "check"))
+  Mesh 2;
+EndIf
+
+For i In {0:Nx_DOM-1}
+For j In {0:Ny_DOM-1}
+  idom = j*Nx_DOM + i;
+  Delete Physicals;
+
+  Physical Point(TAG_CRN~{idom}~{0}) = p~{i}~{j};
+  Physical Point(TAG_CRN~{idom}~{1}) = p~{i+1}~{j};
+  Physical Point(TAG_CRN~{idom}~{2}) = p~{i+1}~{j+1};
+  Physical Point(TAG_CRN~{idom}~{3}) = p~{i}~{j+1};
+  Physical Line(TAG_BND~{idom}~{0}) = lx~{i}~{j};
+  Physical Line(TAG_BND~{idom}~{1}) = ly~{i+1}~{j};
+  Physical Line(TAG_BND~{idom}~{2}) = lx~{i}~{j+1};
+  Physical Line(TAG_BND~{idom}~{3}) = ly~{i}~{j};
+  Physical Surface(TAG_DOM~{idom}) = s~{i}~{j};
+
+  If((i == 0) && (j == 0) && (FLAG_PBM == PBM_SCATT_RECT))
+    Physical Line(TAG_SCA) = {-c1, -c2, -c3, -c4};
+  EndIf
+  If(FLAG_PBM == PBM_MARMOUSI)
+  If((i == I_SOU) && (j == J_SOU))
+    Physical Point(TAG_SOU) = p0;
+  EndIf
+  EndIf
+
+  Printf("Saving subdomain %g...", idom);
+  Save StrCat(MSH_NAME, Sprintf("_%g.msh", idom));
+  Printf("Done.");
+EndFor
+EndFor
+
+Delete Physicals;
+
+If(FLAG_PBM == PBM_SCATT_RECT)
+  Physical Line(TAG_SCA) = {-c1, -c2, -c3, -c4};
+EndIf
+If(FLAG_PBM == PBM_MARMOUSI)
+  Physical Point(TAG_SOU) = p0;
+EndIf
+
+For i In {0:Nx_DOM-1}
+For j In {0:Ny_DOM-1}
+  idom = j*Nx_DOM + i;
+  Physical Surface(TAG_DOM~{idom}) = s~{i}~{j};
+
+  If(j == 0)
+    Physical Line(TAG_BND~{idom}~{0}) = lx~{i}~{j};
+  EndIf
+  If(i == Nx_DOM-1)
+    Physical Line(TAG_BND~{idom}~{1}) = ly~{i+1}~{j};
+  EndIf
+  If(j == Ny_DOM-1)
+    Physical Line(TAG_BND~{idom}~{2}) = lx~{i}~{j+1};
+  EndIf
+  If(i == 0)
+    Physical Line(TAG_BND~{idom}~{3}) = ly~{i}~{j};
+  EndIf
+
+EndFor
+EndFor
+
+Physical Point(TAG_CRN~{          0*Nx_DOM +          0 }~{0}) = p~{0     }~{0     };
+Physical Point(TAG_CRN~{          0*Nx_DOM + (Nx_DOM-1) }~{1}) = p~{Nx_DOM}~{0     };
+Physical Point(TAG_CRN~{ (Ny_DOM-1)*Nx_DOM + (Nx_DOM-1) }~{2}) = p~{Nx_DOM}~{Ny_DOM};
+Physical Point(TAG_CRN~{ (Ny_DOM-1)*Nx_DOM +          0 }~{3}) = p~{0     }~{Ny_DOM};
+
+Printf("Saving domain ...");
+Save StrCat(MSH_NAME, Sprintf(".msh"));
+Printf("Done.");

HelmholtzDDMwithCrosspoints/rectangularDom.pro

+Include "rectangularDom.dat";
+
+// =================================================================================================
+// FUNCTIONS
+// =================================================================================================
+
+Function{
+  If (FLAG_PBM == PBM_SCATT_RECT)
+    f_inc[] = Complex[Cos[WAVENUMBER*X[]], Sin[WAVENUMBER*X[]]];
+    f_ref[] = AcousticFieldSoftCylinder[XYZ[]]{WAVENUMBER, R_SCA};
+    f_dir[] = -f_inc[];
+  EndIf
+  If (FLAG_PBM == PBM_MARMOUSI)
+    f_inc[] = 0;
+    f_ref[] = 0;
+    velocityField[] = InterpolationBilinear[ $1, $2 ]{ ListFromFile["rectangularDomMarmousi.dat"] };
+    cData[] = velocityField[ X[], Y[] ];
+    f_dir[] = -f_inc[];
+    f_sou = 1;
+    For iDom In {0:N_DOM-1}
+      f_sou~{iDom} = 1;
+    EndFor
+  EndIf
+}
+
+// =================================================================================================
+// LISTS ABC/TBC
+// =================================================================================================
+
+// Is there an ABC for main domain? (1=yes, 0=no)
+isEdgeABC~{0} = 1;
+isEdgeABC~{1} = 1;
+isEdgeABC~{2} = 1;
+isEdgeABC~{3} = 1;
+isEdgeRad~{0} = 0;
+isEdgeRad~{1} = 0;
+isEdgeRad~{2} = 0;
+isEdgeRad~{3} = 0;
+
+ListOfSubdom = {};
+For xdom In {0:Nx_DOM-1}
+For ydom In {0:Ny_DOM-1}
+  iDom = ydom*Nx_DOM + xdom;
+
+If (iDom % MPI_Size == MPI_Rank)
+  ListOfSubdom += iDom;
+EndIf
+
+  iDomNeigh~{iDom}~{0} = (ydom-1)*Nx_DOM + xdom;  // Down
+  iDomNeigh~{iDom}~{1} = ydom*Nx_DOM + (xdom+1);  // Right
+  iDomNeigh~{iDom}~{2} = (ydom+1)*Nx_DOM + xdom;  // Top
+  iDomNeigh~{iDom}~{3} = ydom*Nx_DOM + (xdom-1);  // Left
+  iEdgeNeigh~{iDom}~{0} = 2;
+  iEdgeNeigh~{iDom}~{1} = 3;
+  iEdgeNeigh~{iDom}~{2} = 0;
+  iEdgeNeigh~{iDom}~{3} = 1;
+
+  // Is there an ABC or TBC? (1=yes, 0=no)
+  For iEdge In {0:3}
+    isEdgeABC~{iDom}~{iEdge} = 0;
+    isEdgeTBC~{iDom}~{iEdge} = 1;
+    isEdgeRad~{iDom}~{iEdge} = 0;
+  EndFor
+  If(ydom == 0)
+    isEdgeABC~{iDom}~{0} = 1;
+    isEdgeTBC~{iDom}~{0} = 0;
+  EndIf
+  If(xdom == Nx_DOM-1)
+    isEdgeABC~{iDom}~{1} = 1;
+    isEdgeTBC~{iDom}~{1} = 0;
+  EndIf
+  If(ydom == Ny_DOM-1)
+    isEdgeABC~{iDom}~{2} = 1;
+    isEdgeTBC~{iDom}~{2} = 0;
+  EndIf
+  If(xdom == 0)
+    isEdgeABC~{iDom}~{3} = 1;
+    isEdgeTBC~{iDom}~{3} = 0;
+  EndIf
+
+  // List of Edges for Absorbing/Transmission Boundary Condition
+  ListOfEdgesWithABC~{iDom} = {};
+  ListOfEdgesWithTBC~{iDom} = {};
+  For iEdge In {0:3}
+    If(isEdgeABC~{iDom}~{iEdge} == 1)
+      ListOfEdgesWithABC~{iDom} += iEdge;
+    EndIf
+    If(isEdgeTBC~{iDom}~{iEdge} == 1)
+      ListOfEdgesWithTBC~{iDom} += iEdge;
+    EndIf
+  EndFor
+
+EndFor
+EndFor
+
+// =================================================================================================
+// GROUPS
+// =================================================================================================
+
+Group {
+
+// === REFERENCE CASE ===
+
+  // DOMAIN
+  Omega = Region[{}];
+
+  // SOURCE
+  GammaPoint = Region[{}];
+
+  // BOUNDARY (Dirichlet BC)
+  GammaD0 = Region[{}];
+  GammaN0 = Region[{}];
+  GammaD = Region[{}];
+  GammaN = Region[{}];
+
+  // BOUNDARY (Absorbing BC)
+  For iEdge In {0:3}
+    SigmaMain~{iEdge} = Region[{}];
+    SigmaMainSide~{iEdge}~{0} = Region[{}];
+    SigmaMainSide~{iEdge}~{1} = Region[{}];
+  EndFor
+  SigmaMain = Region[{}];
+
+  // CORNER (Absorbing BC)
+  CornerMain~{0} = Region[{TAG_CRN~{          0*Nx_DOM +          0 }~{0}}];
+  CornerMain~{1} = Region[{TAG_CRN~{          0*Nx_DOM + (Nx_DOM-1) }~{1}}];
+  CornerMain~{2} = Region[{TAG_CRN~{ (Ny_DOM-1)*Nx_DOM + (Nx_DOM-1) }~{2}}];
+  CornerMain~{3} = Region[{TAG_CRN~{ (Ny_DOM-1)*Nx_DOM +          0 }~{3}}];
+  CornerMain = Region[{ CornerMain~{0}, CornerMain~{1}, CornerMain~{2}, CornerMain~{3}}];
+
+  If(FLAG_PBM == PBM_SCATT_RECT)
+    GammaD += Region[{TAG_SCA}];
+  EndIf
+
+  If(FLAG_PBM == PBM_MARMOUSI)
+    GammaPoint += Region[{TAG_SOU}];
+  EndIf
+
+For xdom In {0:Nx_DOM-1}
+For ydom In {0:Ny_DOM-1}
+  iDom = ydom*Nx_DOM + xdom;
+
+  // DOMAIN
+  Omega += Region[TAG_DOM~{iDom}];
+
+  // BOUNDARY (Absorbing BC)
+  For iEdge In {0:3}
+  If(ydom == 0)
+    SigmaMain~{0} += Region[{TAG_BND~{iDom}~{0}}];
+    SigmaMain += Region[{TAG_BND~{iDom}~{0}}];
+  EndIf
+  If(xdom == Nx_DOM-1)
+    SigmaMain~{1} += Region[{TAG_BND~{iDom}~{1}}];
+    SigmaMain += Region[{TAG_BND~{iDom}~{1}}];
+  EndIf
+  If(ydom == Ny_DOM-1)
+    SigmaMain~{2} += Region[{TAG_BND~{iDom}~{2}}];
+    SigmaMain += Region[{TAG_BND~{iDom}~{2}}];
+  EndIf
+  If(xdom == 0)
+    SigmaMain~{3} += Region[{TAG_BND~{iDom}~{3}}];
+    SigmaMain += Region[{TAG_BND~{iDom}~{3}}];
+  EndIf
+  EndFor
+EndFor
+EndFor
+
+  For iEdge In {0:3}
+    iCornerSide1 = iEdge;
+    iCornerSide2 = (iEdge+1)%4;
+    SigmaMainSide~{iEdge}~{0} = Region[{CornerMain~{iCornerSide1}}];
+    SigmaMainSide~{iEdge}~{1} = Region[{CornerMain~{iCornerSide2}}];
+  EndFor
+
+// === DDM CASE ===
+
+For i0 In {0:#ListOfSubdom()-1}
+iDom = ListOfSubdom(i0);
+
+  // DOMAIN
+  Omega~{iDom} = Region[TAG_DOM~{iDom}];
+
+  // SOURCE
+  GammaPoint~{iDom} = Region[{}];
+  If(FLAG_PBM == PBM_MARMOUSI)
+  If(iDom == Idom_SOU)
+    GammaPoint~{iDom} += Region[{TAG_SOU}];
+  EndIf
+  EndIf
+
+  // BOUNDARY (Dirichlet/Neumann BC)
+  GammaD0~{iDom} = Region[{}];
+  GammaN0~{iDom} = Region[{}];
+  GammaD~{iDom} = Region[{}];
+  GammaN~{iDom} = Region[{}];
+  If((FLAG_PBM == PBM_SCATT_RECT) && (iDom == 0))
+    GammaD~{iDom} = Region[{TAG_SCA}];
+  EndIf
+  TrOmegaGammaD~{iDom} = ElementsOf[ Omega~{iDom}, OnOneSideOf GammaD~{iDom} ];
+
+  // BOUNDARY (Absorbing BC + Transmission BC)
+  For iEdge In {0:3}
+    iCornerSide1 = iEdge;
+    iCornerSide2 = (iEdge+1)%4;
+    Sigma~{iDom}~{iEdge} = Region[{TAG_BND~{iDom}~{iEdge}}];
+    SigmaSide~{iDom}~{iEdge}~{0} = Region[{TAG_CRN~{iDom}~{iCornerSide1}}];
+    SigmaSide~{iDom}~{iEdge}~{1} = Region[{TAG_CRN~{iDom}~{iCornerSide2}}];
+    GammaD~{iDom}~{iEdge}  = Region[{}];
+    GammaD0~{iDom}~{iEdge} = Region[{}];
+    GammaN~{iDom}~{iEdge}  = Region[{}];
+    GammaN0~{iDom}~{iEdge} = Region[{}];
+  EndFor
+  Sigma~{iDom} = Region[{ Sigma~{iDom}~{0}, Sigma~{iDom}~{1}, Sigma~{iDom}~{2}, Sigma~{iDom}~{3} }];
+
+  // CORNER (Absorbing BC + Transmission BC)
+  For iCorner In {0:3}
+    Corner~{iDom}~{iCorner} = Region[{TAG_CRN~{iDom}~{iCorner}}];
+  EndFor
+  Corner~{iDom} = Region[{ Corner~{iDom}~{0}, Corner~{iDom}~{1}, Corner~{iDom}~{2}, Corner~{iDom}~{3} }];
+
+EndFor
+}