Merge branch 'master' of https://gitlab.onelab.info/doc/models

HelmholtzDDMwithCrosspoints/README.txt

+A non-overlapping HABC-based domain decomposition method (DDM) with
+cross-point treatment for 2D Helmholtz problems.
+
+Models developed by Axel Modave, based on GetDDM codes.
+
+
+Quick start
+-----------
+
+Open `main.pro' with Gmsh.
+
+
+Additional info
+---------------
+
+Used in:
+
+  A. Modave, A. Royer, X. Antoine, X. Geuzaine. An optimized Schwarz domain
+  decomposition method with cross-point treatment for time-harmonic acoustic
+  scattering.  Preprint: https://hal.archives-ouvertes.fr/hal-02432422

HelmholtzDDMwithCrosspoints/domCircle.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/domCircle.geo

+Include "domCircle.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/domCircle.pro

+Include "domCircle.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/domRectangle.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/domRectangle.geo

+Include "domRectangle.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/domRectangle.pro

+Include "domRectangle.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["domRectangleMarmousi.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
+}

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 = "domCircle.geo";
+  LinkPro = "domCircle.pro";
+EndIf
+If ((FLAG_PBM == PBM_SCATT_RECT) ||  (FLAG_PBM == PBM_MARMOUSI))
+  LinkGeo = "domRectangle.geo";
+  LinkPro = "domRectangle.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];

HelmholtzHABCwithCorners/README.txt

+High-order absorbing boundary conditions (HABCs) with corner and edge treatment
+for 2D and 3D Helmholtz problems.
+
+Models developed by Axel Modave.
+
+
+Quick start
+-----------
+
+Open `main.pro' with Gmsh.
+
+
+Additional info
+---------------
+
+Used in:
+
+  A. Modave, X. Geuzaine, X. Antoine (2020). Corner treatments for high-order
+  absorbing boundary conditions in high-frequency acoustic scattering problems.
+  J. Comput. Phys., 401, 109029.
+  Preprint: https://hal.archives-ouvertes.fr/hal-01925160

HelmholtzHABCwithCorners/domCube.dat

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

HelmholtzHABCwithCorners/domCube.geo

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

HelmholtzHABCwithCorners/domDisk.dat

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

HelmholtzHABCwithCorners/domDisk.geo

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

HelmholtzHABCwithCorners/domDisk.pro

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