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

AcousticScattering/.gitignore

+*.pos
+*.msh
+*.pre
+*.db

DDM-Corner-Helmholtz-2D/.gitignore

+*.pos
+*.msh
+*.pre
+*.db
+out_ABC/**
+out_TC/**
\ No newline at end of file

DDM-Corner-Helmholtz-2D/ABC_norm.py

+import sys
+import os, os.path
+import numpy as np
+import subprocess
+import matplotlib.pyplot as plt
+
+def sumNorm(foldername, rootname):
+  # Find number of files
+  nfile = 0
+  for name in os.listdir(foldername):
+    if os.path.isfile(os.path.join(foldername, name)):
+      if rootname in name:
+        nfile += 1
+  print("Number of files : " + str(nfile))
+
+  normL2 =[]
+  for n in range(nfile):
+    norm = 0.
+    f = open(os.path.join(foldername, rootname +str(n)), "r")
+    for line in f:
+      if(len(line) != 1):
+        values = line.split()
+        norm += float(values[1])
+    f.close()
+    normL2.append(np.sqrt(norm))
+
+  print(normL2)
+  return(normL2)
+
+foldername = "out_ABC/"
+rootname = "ABC_DDM_normg.pos"
+NDom = 6
+MaxIt = 100
+
+subprocess.run(["rm", "-rf", foldername])
+subprocess.run(["gmsh", "ABC_main.geo", "-2"])
+subprocess.run(["getdp", "main.pro", "-solve", "DDM_1",
+                        "-v", "3",
+                        "-setstring", "PROBLEM", "Absorbing BC",
+                        "-setstring", "SOLVER", "jacobi",
+                        "-setstring", "ABC_DIRREL", foldername,
+                        "-setnumber", "MAXIT", str(MaxIt),
+                        "-setnumber", "Ndom", str(NDom)])
+
+errorL2 = sumNorm(foldername, rootname)
+
+plt.plot(errorL2)
+plt.show()
\ No newline at end of file

DDM-Corner-Helmholtz-2D/README.md

+# ABC and Transmission CondtionCorner
+
+## Requierement
+
+Two open-source software are requiered to be installed:
+
+- [GMSH](https://gmsh.info), *latest automatic Gmsh snapshot* (Version > 4.5.6 (Not the 4.5.6))
+- [GetDP](http://getdp.info/) : Version >=  3.3.0
+
+GMSH is a mesh generator with post-processing facilities while GetDP is a finite element solver.
+
+## How to use
+
+Launch the script `main.pro` with GMSH, either through the GUI (Graphical User Interface) or from CLI (Command Line Interface):
+```bash
+gmsh main.pro
+```
+
+On the left-sided menu, choose the `Problem` ("Absorbing BC" or "Transmission Condition"), set the `Resolution` and the different parameters in the `Input` sub-menu and finally click on `Run`. A consol log can be show by clicking on the bottom line of the windows on the left part.
+
+> See the wiki for an animated summary!
+
+## Parameters and Options
+
+### ABC
+
+The script produces three functions:
+- `u`: Scattered field (not the total field!)
+- `uNorm`: Absolute value of the scattered field
+- `err`: Error (difference) between `u` and an analytic solution computed thanks to Mie Serie decomposition
+
+
+![Comparaison of result with `DDM-1` algorithm, with 2 conditions at corners: homog. neumann and our corner correction.](img/ABC.png)
+
+Here are the main options that can be tuned by the user:
+
+- `ABC`: Order of the condition (0 or 2) and the type of corner condition (Homogeneous Neumann or with our "corner correction")
+- `DDM`: Number of subdomain: 2,3 or 6 (if launched with a DD-Algorithm). The monodomain problem can  be launched independently of the geometry (see section `GetDP`). The `beta` parameter cannot be changed.
+- `Geometry`: the radii of the ABC and obstacle and the position of the obstacle.
+- `GetDP`: The `Resolution` value:
+  - `MonoDomain`: Direct solver / Mono domain problem (every sub-domains are merged)
+  - `DDM-1`: First algorithm presented in the paper: auxiliary function `phi` is global
+  - `DDM-2`: Second algorithm where `phi` is local to each subdomain
+- `GMSH`: nothing
+- `Input` : `Incident angle` of the plane wave and `Wavenumber`
+- `IterativeSolver`
+  - `Solver`: `Jacobi` ("Parallel Schwarz") or `gmres`. The `print` is only for debugging purpose and might not work on your configuration.
+  - `Tolerance`, `Max it` and `Restart` (for GMRES) are classical parameter
+- `Mesh`: `NLambda` is the number of discretization points per wavelength
+- `Output`: `Output Directory` of the results
+
+## Transmission Condition
+
+The script returns the total field `u` in the waveguide.
+
+![The transmission condition problem: a waveguide divided in 2 subdomains with a possibly broken line as a border.](img/TC.png)
+
+Here are the available options:
+
+- `GetDP`: The `Resolution` value:
+  - `MonoDomain`: Direct solver / Mono domain problem (every sub-domains are merged)
+  - `DDM`: Domain Decomposition Algorithm
+- `GMSH`: nothing
+- `DDM`: 
+  - `Order 2` (continuous auxiliary function), `Homog. Neumann` or (our) `Corner Correction`
+  - `Corner Condition`: `Dirichlet` (continuous auxiliary function), `Homog. Neumann` or (our) `Corner Correction`
+- `Geometry`: 
+  - `X-width` and `Y-width`: resp. X-length and Y-length of the waveguide
+  - `Type of border line`: `Broken Line` or `Straight Line`
+  - `X-coord of the pick point`: Move the middle (peak) point on the x-line
+  - `Y-coord of bottom point`: Move the bottom point on the y-line. The top point is moved symmetrically
+- `Boundary Conditions`:
+  - `Incoming (left)`: the boundary condition of the left side of the square, where the incoming wave is sent. It can be either `Fourier` (dn u + iwu) or `Dirichlet`.
+  - `Outgoing (right, =0)`: either `Fourier` (dn u + iwu), `Neumann` or `Dirichlet`. The condition is homogeneous (=0) and set on the right side of the square.
+  - `Top (=0)`: either `Neumann` or `Dirichlet`, homogeneous in both case (=0), for the top side of the square.
+  - `Bottom (=0)`: either `Neumann` or `Dirichlet`, homogeneous in both case (=0), for the bottom side of the square.
+- `Input`: 
+  - `wavenumber`
+  - `Type of incident wave`: `Plane wave` (exp^{i*k*(alpha*x)}, alpha = `Incident angle`) or Fourier mode (with m = `Mode number`)
+- `IterativeSolver`
+  - `Solver`: `Jacobi` ("Parallel Schwarz") or `gmres`. The `print` is only for debugging purpose and might not work on your configuration.
+  - `Tolerance`, `Max it` and `Restart` (GMRES only) are classical parameter
+- `Mesh`: `NLambda` is the number of discretization points per wavelength. The quantity `h` is the diameter of an element.
+- `Output`: 
+  - `Prefix for filename`: prefix applied to every saved file
+  - `Output Directory` of the results
+  - `Print every phi`: print on disk every auxiliary functions
+
+# Authors
+
+Bruno Després, Anouk Nicolopoulos, Bertrand Thierry
+# License
+
+GNU General Public License v3.0 or later
+
+See COPYING to see the full text.
\ No newline at end of file

DDM-Corner-Helmholtz-2D/main.data

+DefineConstant[
+  PROBLEM = {"Absorbing BC", Choices {"Absorbing BC", "Transmission Condition"}, Visible 1, Name "GetDP/0Problem", GmshOption "Reset", Highlight "Blue"}
+];
+
+
+// GetDP
+DefineConstant[
+  // default getdp parameters for onelab
+  R_ = {"MonoDomain", Name "GetDP/1ResolutionChoices", Visible 1, Highlight "Blue"},
+  C_ = {"-solve -v 3 -bin -ksp_monitor", Name "GetDP/9ComputeCommand", Visible 0},
+  P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}
+];
\ No newline at end of file

DDM-Corner-Helmholtz-2D/main.geo

+WorkingDir = CurrentDir;
+Include "main.data";
+
+If(!StrCmp(PROBLEM, "Absorbing BC"))
+  Include "src/ABC_main.data";
+  Include "src/ABC_main.geo";
+Else
+  Include "src/TC_main.data";
+  Include "src/TC_main.geo";
+EndIf

DDM-Corner-Helmholtz-2D/main.pro

+WorkingDir = CurrentDir;
+Include "main.data";
+
+Include "src/Common.pro";
+
+If(!StrCmp(PROBLEM,"Absorbing BC"))
+  Include "src/ABC_main.data";
+  Include "src/ABC_main.pro";
+Else
+  Include "src/TC_main.data";
+  Include "src/TC_main.pro";
+EndIf
+
+

DDM-Corner-Helmholtz-2D/notes.md

+# Notations
+
+## GetDP Group (`hexa.pro`)
+
+Les indices `i` et `j` sont réservés aux sous domaines tandis que `k` et `l` le sont pour les segments du bord extérieur.
+
+### Numbers
+
+| Nom| Fonction|
+| ----- | ----- |
+|`NGammaExt` |Number of segments of exterior boundary $\Gamma$ |
+|`NGammaExt~{i}` |Number of segments of $\Gamma\cap\partial\Omega_i$ |
+|`NDom` |  Number of subdomains|
+
+### Indices
+
+| Nom| Fonction|
+| ----- | ----- |
+|`Omega~{i}` | Domaine $\Omega_i$|
+|`GammaExt` | Bord $\Gamma$ en entier|
+|`GammaExt~{k}` | Segment $\Gamma_k$|
+| `GammaExt~{i}~{k}`| Segment $\Gamma_k^i$|
+|`GammaInt` | Bord intérieur (obstacle) en entier|
+|`GammaInt~{i}` | `GammaInt`$\cap\partial\Omega_i$|
+| `Apt~{i}~{j}~{k}~{l}` | Point $A_{k\ell}^{ij}$ |
+| `dGammaExt~{i}~{k}~{0}` | Same as one of the `Apt~{i}~{j}~{k}~{l}` but listed differently |
+| `Cset~{i}~{k}` | End-points of $\Gamma_k^i$|
+| `dGammaExt~{i}~{k}` | = `Cset~{i}~{k}`|
+| `dOmegaInt~{i}` | $=$`GammaInt~{i}`$=\partial\Omega_i \cap$`GammaInt`|
+| `dOmegaExt~{i}` | $=\partial\Omega_i \cap$`GammaExt`$=\cup_k\Gamma_k^i$|
+| `Sigma~{i}~{j}` | $=\Sigma_{ij} = \partial\Omega_i\cap\partial\Omega_j$|
+| `Sigma~{i}` | $=\cap_j$ `Sigma~{i}~{j}`|
+| `dSigmaExt~{i}~{j}~{n}` <br> with `n = 0,1,2...` | $n^{th}$ (End-) Point of $\Sigma_{ij}$ such that it also belons to $\Gamma$ (could be more than 1 end-point !)|
+| `dSigmaInt~{i}~{j}~{n}` <br> with `n = 0,1,2...` | Same but switch the role of `GammaExt`($=\Gamma$) with `GammaInt`|
+| `dSigmaExt~{i}~{j}` | $=\partial\Sigma_{ij}\cap\Gamma= \cup_n$`dSigmaExt~{i}~{j}~{n}`|
+| `dSigmaInt~{i}~{j}` |  $=\partial\Sigma_{ij}\cap\Gamma^{int}= \cup_n$`dSigmaInt~{i}~{j}~{n}`|
+| `dSigma~{i}~{j}` |  $=\partial\Sigma_{ij}=$`dSigmaExt~{i}~{j}`$\cup$`dSigmaInt~{i}~{j}`|
+
+### Lists
+
+> A list in GetDP is always monodimensional. When a list is said to be of dimension 2, explanation are provided to explain how to access the data.
+
+| Nom| Dim | Fonction|
+| ----- | ----- |----- |
+|`D()` | 1 |Indices $i$ of subdomains $\Omega_{\texttt{D(index)}}$ with `index=0,1,2...,Ndom-1`. The index `i` of the subdomain is obtained through: `i= D(index)`. When subdomains are numbered from 0 to `Ndom`-1, then `D()` is simply equal to `[0,1,2,...Ndom-1]`.|
+|`myD()` | 1 |Indices `index=0,1,2,...` in `D()` such that $\Omega_{\texttt{D(index)}}$ is  a subdmain the current MPI process has to treat (sequential case: `myD() = 0,1,2,...,Ndom-1`)|
+| `D~{i}`| 1 |Indices $j$ of neighbor : *i.e.* such that $\partial\Omega_j\cap\partial\Omega_i\neq\emptyset$. Note that `i` and `j` are the "real" subdomain numbers |
+| `ListGammaExt~{i}()` | 1 |Indices $k$ (global) of $\Gamma_k^i$ ($\neq \emptyset$ by definition) |
+| `ListdGammaExt~{i}~{k}()` | 2 | Couple of indices $(j, \ell)$ of the $n^{th}$ end-point of $\Gamma_k^i$ (*i.e.* the $j$ and $\ell$ of $A_{k\ell}^{ij}$. Store contiguously : `[j,l, j,l, j,l,...]` |
+| `ListdSigmaExt~{i}~{j}()` | 2 | Couple of indices $(k, \ell)$ of the $n^{th}$ end-point of $\Sigma_{ij}$ such that it also belongs to $\Gamma$ (*i.e.* the $k$ and $\ell$ of $A_{k\ell}^{ij}$. Store contiguously : `[j,l, j,l, j,l,...]` (same role as `ListdGammaExt~{i}~{k}()`) | 
+

DDM-Corner-Helmholtz-2D/src/ABC_DDM.pro

+// DDM-algorithms
+// Common quantities and operations are defined here
+// DDM algorithm are included below 
+
+// Numbering of the unknowns of the DD algorithms
+Include "ABC_Decomposition.pro";
+
+Group {
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    TrdOmegaInt~{i} = ElementsOf[Omega~{i}, OnOneSideOf dOmegaInt~{i}];
+  EndFor
+}
+
+Constraint {
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    {Name Dirichlet~{i}; Type Assign; Case{{Region dOmegaInt~{i}; Value $PhysicalSource?-uInc[]:0.; } } }
+  EndFor
+}
+
+FunctionSpace{
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    {Name H_grad~{i}; Type Form0;
+      BasisFunction{
+        {Name u~{i}; NameOfCoef ui; Function BF_Node;
+      	Support Region[{Omega~{i}, dOmegaExt~{i}, dOmegaInt~{i}, Sigma~{i}}]; Entity NodesOf[All];}
+        }
+      Constraint{
+        {NameOfCoef ui; EntityType NodesOf; NameOfConstraint Dirichlet~{i};}
+      }
+    }
+
+    For jj In {0:#myD~{i}()-1}
+      j = myD~{i}(jj);
+      {Name H_g~{i}~{j}; Type Form0;
+        BasisFunction{
+        {Name g~{i}~{j}; NameOfCoef gj; Function BF_Node;
+	        Support Region[{Sigma~{i}~{j}}];Entity NodesOf[Sigma~{i}~{j}];}
+        }
+      }
+    EndFor
+  EndFor
+}
+
+// ============ DDM Algorithms
+Include "ABC_DDM_1.pro";
+Include "ABC_DDM_2.pro";
+// ===============
+
+// Common Post-Processings to every DDM algorithms
+// ===============================================
+PostProcessing{
+  // Transmitted Data
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    For jj In {0:#myD~{i}()-1}
+      j = myD~{i}(jj);
+      {Name G~{j}~{i}; NameOfFormulation G~{j}~{i};
+        Quantity {
+          {Name g~{j}~{i}; Value {Local { [{g~{j}~{i}}] ; In Sigma~{i}~{j}; Jacobian JSur; }}}
+        }
+      }
+      {Name NormL2G~{j}~{i}; NameOfFormulation G~{j}~{i};
+        Quantity {
+          {Name normL2g~{j}~{i}; Value { Integral { [ Abs[g~{j}~{i}[] - {g~{j}~{i}} - b~{j}~{i}[]]^2 ]; In Sigma~{i}~{j}; Jacobian JSur; Integration I1; }}}
+        }
+      }
+    EndFor
+  EndFor
+  // Local solutions (whatever the algorithm)
+  {Name DDM; NameOfFormulation DDM_1;
+    Quantity {
+      {Name u_ddm; Value {
+    	  For ii In {0:#myD()-1}
+	        i = myD(ii);
+	        Local { [{u~{i}}] ; In Omega~{i}; Jacobian JVol; }
+	      EndFor
+    	  }
+      }
+      {Name u_ddm_abs; Value {
+	      For ii In {0:#myD()-1}
+	        i = myD(ii);
+	        Local { [Norm[{u~{i}}]] ; In Omega~{i}; Jacobian JVol; }
+	      EndFor
+      	}
+      }
+      //Error between solution and exact solution (exact in free space!)
+      {Name uerr_ddm; Value {
+	      For ii In {0:#myD()-1}
+	        i = myD(ii);
+	        Local { [Norm[{u~{i}} - uex[]]] ; In Omega~{i}; Jacobian JVol; }
+	      EndFor
+    	  }
+      }
+    }
+  }
+}
+
+// Common Post-Operations to every DDM-algorithms
+// ===============================================
+PostOperation{
+  // Compute the transmitted data + Store it in memory
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    For jj In {0:#myD~{i}()-1}
+      j = myD~{i}(jj);
+      {Name G~{j}~{i}; NameOfPostProcessing G~{j}~{i};
+      	Operation { Print [g~{j}~{i}, OnElementsOf Sigma~{i}~{j}, StoreInField tag_g~{j}~{i}];}
+      }
+      //Store RHS to compute L2 error (Jacobi)
+      {Name B~{j}~{i}; NameOfPostProcessing G~{j}~{i};
+      	Operation { Print [g~{j}~{i}, OnElementsOf Sigma~{i}~{j}, StoreInField tag_b~{j}~{i}];}
+      }
+      {Name NormL2G~{j}~{i}; NameOfPostProcessing NormL2G~{j}~{i};
+        Operation{
+          Print [normL2g~{j}~{i}[Sigma~{i}~{j}], OnRegion Sigma~{i}~{j}, Format Table, AppendExpressionToFileName $ItCount, File >> StrCat[ABC_DIR,"ABC_DDM_normg.pos"]];
+        }
+      }
+    EndFor
+  EndFor
+  // Print the solution on HDD + error wrt exact solution
+  {Name DDM; NameOfPostProcessing DDM;
+    Operation {
+      Print [u_ddm, OnElementsOf Omega, File StrCat[ABC_DIR, StrCat[ABC_PREFIX, "DDM_Scattered_Field.pos"]]];
+      Print [u_ddm_abs, OnElementsOf Omega, File StrCat[ABC_DIR, StrCat[ABC_PREFIX, "DDM_Scattered_Field_(abs._value).pos"]]];
+      Print [uerr_ddm, OnElementsOf Omega, File StrCat[ABC_DIR, StrCat[ABC_PREFIX, "DDM_Error_numeric_vs_analytic.pos"]]];
+    }
+  }
+}

DDM-Corner-Helmholtz-2D/src/ABC_DDM_1.pro

+// First algorithm
+
+Formulation {
+  { Name DDM_1; Type FemEquation;
+    Quantity{
+			For ii In {0:#myD()-1}
+				i = myD(ii);
+				{ Name u~{i} ; Type Local; NameOfSpace H_grad~{i};}
+			EndFor
+			For k In {0:#ListGammaExt()-1}
+				{ Name phi~{k} ; Type Local; NameOfSpace H_phi~{k};}
+			EndFor
+		}
+    Equation{
+    	For ii In {0:#myD()-1}
+				i = myD(ii);
+				// PDE on u_i
+				//Helmholtz equation
+				Galerkin{[Dof{Grad u~{i}}, {Grad u~{i}}];
+					In Omega~{i}; Jacobian JVol; Integration I1;}
+				Galerkin{[-kap^2*Dof{u~{i}}, {u~{i}}];
+					In Omega~{i}; Jacobian JVol; Integration I1;}
+				// dni ui -i*kap ui = gij on Sigma ij
+				For jj In {0:#myD~{i}()-1}
+					j = myD~{i}(jj);
+					Galerkin{[-I[]*kap*Dof{u~{i}}, {u~{i}}];
+							In Sigma~{i}~{j}; Jacobian JSur; Integration I1;}
+					Galerkin{[$TCSource?-g~{i}~{j}[]:0., {u~{i}}];
+						In Sigma~{i}~{j}; Jacobian JSur; Integration I1;}
+				EndFor
+				If(ABC_ORDER==0)
+					Galerkin{[-I[]*kap*Dof{u~{i}}, {u~{i}}];
+						In dOmegaExt~{i}; Jacobian JSur; Integration I1;}
+				Else
+					// dn u_i = i*omega*phi_k on Gamma_i^k
+					For kk In {0:#ListGammaExt~{i}()-1}
+						k=ListGammaExt~{i}(kk);
+						Galerkin{[-I[]*kap*Dof{phi~{k}}, {u~{i}}];
+							In GammaExt~{i}~{k}; Jacobian JSur; Integration I1;}
+					EndFor
+				EndIf
+			EndFor
+
+			// Operator T
+			If(ABC_ORDER == 2)
+				For ii In {0:#myD()-1}
+					i = myD(ii);
+					For kk In {0:#ListGammaExt~{i}()-1}
+						k = ListGammaExt~{i}(kk);
+						Galerkin{[-Dof{u~{i}}, {phi~{k}}];
+							In GammaExt~{i}~{k}; Jacobian JSur; Integration I1;}
+					EndFor
+				EndFor
+				For kk In {0:#ListGammaExt()-1}
+					k = ListGammaExt(kk);
+					Galerkin{[Dof{phi~{k}}, {phi~{k}}];
+						In GammaExt~{k}; Jacobian JSur; Integration I1;}
+					Galerkin{[1/(2*kap^2)*Dof{d phi~{k}}, {d phi~{k}}];
+						In GammaExt~{k}; Jacobian JSur; Integration I1;}
+					// Corner correction
+					If(CORNER_CONDITION == 2)
+						For index In {0:#ListdGammaExt~{k}()-1}
+							// Loop on end-point
+							l=ListdGammaExt~{k}(index); // Index of ajacent segment
+							Galerkin{[-I[]/(4*kap)*(Cos[theta]/Cos[theta/2]+Cos[theta/2])*Dof{phi~{k}}, {phi~{k}}];
+								In Apt~{k}~{l}; Jacobian JLin; Integration I1;}
+							Galerkin{[-I[]/(4*kap)*(-Cos[theta]/Cos[theta/2]+Cos[theta/2])*Dof{phi~{l}}, {phi~{k}}];
+								In Apt~{k}~{l}; Jacobian JLin; Integration I1;}
+						EndFor
+					EndIf
+				EndFor
+			EndIf
+    }
+  }
+
+  // Equation to compute the transmitted data g_{ji} that is sent to Omega_j
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    For jj In {0:#myD~{i}()-1}
+    	j = myD~{i}(jj);
+      	{ Name G~{j}~{i}; Type FemEquation;
+	        Quantity{
+						{ Name u~{i} ; Type Local; NameOfSpace H_grad~{i};}
+						{ Name g~{j}~{i} ; Type Local; NameOfSpace H_g~{j}~{i};}
+					}
+			Equation{
+				Galerkin{[Dof{g~{j}~{i}},{g~{j}~{i}}];
+					In Sigma~{j}~{i}; Jacobian JSur; Integration I1;}
+				Galerkin{[$TCSource?+g~{i}~{j}[]:0.,{g~{j}~{i}}];
+	    		 In Sigma~{j}~{i}; Jacobian JSur; Integration I1;}
+				Galerkin{[2*I[]*kap*{u~{i}},{g~{j}~{i}}];
+	    		 In Sigma~{j}~{i}; Jacobian JSur; Integration I1;}
+				}
+			}
+    EndFor
+  EndFor
+}
+
+Resolution{
+  {Name DDM_1;
+    System{
+      {Name DDM_1; NameOfFormulation DDM_1; Type Complex; }
+      For ii In {0:#myD()-1}
+				i = myD(ii);
+				For jj In {0:#myD~{i}()-1}
+	  			j = myD~{i}(jj);
+	  			{Name G~{j}~{i}; NameOfFormulation G~{j}~{i}; Type Complex; }
+				EndFor
+      EndFor
+    }
+    Operation{
+			CreateDir[ABC_DIR];
+      // Initialisation
+      Evaluate[$PhysicalSource = 1];
+      Evaluate[$TCSource = 0];
+      UpdateConstraint[DDM_1, dOmegaInt, Assign];
+      Generate[DDM_1];
+      Solve[DDM_1];
+      For ii In {0:#myD()-1}
+				i = myD(ii);
+				For jj In {0:#myD~{i}()-1}
+	  			j = myD~{i}(jj);
+	  			Generate[G~{j}~{i}];
+	  			Solve[G~{j}~{i}];
+					PostOperation[G~{j}~{i}];
+					//Save RHS if jacobi
+					If(StrCmp(SOLVER, "jacobi") == 0)
+						PostOperation[B~{j}~{i}];
+					EndIf
+				EndFor
+      EndFor
+      Evaluate[$PhysicalSource = 0];
+      Evaluate[$TCSource = 1];
+			UpdateConstraint[DDM_1, dOmegaInt, Assign];
+			Evaluate[$ItCount = 0];
+      IterativeLinearSolver["I-A",SOLVER, TOL, MAXIT, RESTART, {ListOfFields_g()}, {}, {}]
+      {
+				// Compute local solutions
+				GenerateRHSGroup[DDM_1, Region[{Sigma, TrdOmegaInt, GammaExt, dGammaExt}]];
+				SolveAgain[DDM_1];
+				For ii In {0:#myD()-1}
+				  i = myD(ii);
+				  For jj In {0:#myD~{i}()-1}
+				    j = myD~{i}(jj);
+				    GenerateRHSGroup[G~{j}~{i}, Region[{Sigma~{i}}]];
+				    SolveAgain[G~{j}~{i}];
+				  EndFor
+				EndFor
+				// Update data to be exchanged
+				For ii In {0:#myD()-1}
+					i = myD(ii);
+					For jj In {0:#myD~{i}()-1}
+						j = myD~{i}(jj);
+						If(StrCmp(SOLVER, "jacobi") == 0)
+							PostOperation[NormL2G~{j}~{i}]; // Compute difference between new and old g_ji
+						EndIf
+						PostOperation[G~{j}~{i}];
+					EndFor
+				EndFor
+
+				Evaluate[$ItCount = $ItCount+1];
+      }
+      //Compute global solution
+      Evaluate[$PhysicalSource = 1];
+      Evaluate[$TCSource = 1];
+      UpdateConstraint[DDM_1, dOmegaInt, Assign];
+      GenerateRHSGroup[DDM_1, Region[{Sigma, TrdOmegaInt, GammaExt, dGammaExt}]];
+      SolveAgain[DDM_1];
+      PostOperation[DDM];
+    }
+  }
+}

DDM-Corner-Helmholtz-2D/src/ABC_DDM_2.pro

+// Algorithm DDM-2
+Function{
+  coef_sigma[] = -(1 + I[]*beta/kap*Cos[theta/2]/Cos[theta]);
+  coef_gamma[] = beta + I[]*kap*Cos[theta/2];
+  If(ABC_ORDER==0)
+    ListOfFields_ddm2() = {ListOfFields_g()};
+  Else
+    ListOfFields_ddm2() = {ListOfFields_g(),ListOfFields_h()};
+  EndIf
+}
+
+FunctionSpace{
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    For kk In {0:#ListGammaExt~{i}()-1}
+      k=ListGammaExt~{i}(kk);
+      {Name H_phi~{i}~{k}; Type Form0;
+        BasisFunction{
+          {Name phi~{i}~{k}; NameOfCoef phiik; Function BF_Node;
+            Support Region[{GammaExt~{i}~{k}, dGammaExt~{i}~{k}}];
+            Entity NodesOf[All];}
+        }
+      }
+      // Quantities transfered at each Apt~{i}~{j}~{k}~{l}
+      For index In {0:#ListdGammaExt~{i}~{k}()-1:2}
+        j=ListdGammaExt~{i}~{k}(index);
+        l=ListdGammaExt~{i}~{k}(index+1);
+        {Name H_h~{i}~{j}~{k}~{l}; Type Form0;
+          BasisFunction{
+            {Name h~{i}~{j}~{k}~{l}; NameOfCoef hj; Function BF_Node; 
+            Support Region[{Apt~{i}~{j}~{k}~{l}}]; Entity NodesOf[All];}
+          }
+        }
+      EndFor
+    EndFor
+  EndFor
+}
+
+Formulation{
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    { Name DDM_2~{i}; Type FemEquation;
+      Quantity{
+      	{Name u~{i} ; Type Local; NameOfSpace H_grad~{i};}
+        For kk In {0:#ListGammaExt~{i}()-1}
+          k=ListGammaExt~{i}(kk);
+          { Name phi~{i}~{k} ; Type Local; NameOfSpace H_phi~{i}~{k};}
+        EndFor
+      }
+      Equation{
+        // PDE on u_i
+				//Helmholtz equation
+        Galerkin{[Dof{Grad u~{i}}, {Grad u~{i}}];
+          In Omega~{i}; Jacobian JVol; Integration I1;}
+        Galerkin{[-kap^2*Dof{u~{i}}, {u~{i}}];
+          In Omega~{i}; Jacobian JVol; Integration I1;}
+        // dni ui -i*kap ui = gij on Sigma ij
+        For jj In {0:#myD~{i}()-1}
+          j = myD~{i}(jj);
+          Galerkin{[-I[]*kap*Dof{u~{i}}, {u~{i}}];
+            In Sigma~{i}~{j}; Jacobian JSur; Integration I1;}
+          Galerkin{[$TCSource?-g~{i}~{j}[]:0., {u~{i}}];
+            In Sigma~{i}~{j}; Jacobian JSur; Integration I1;}
+        EndFor
+        If(ABC_ORDER==0)
+					Galerkin{[-I[]*kap*Dof{u~{i}}, {u~{i}}];
+						In dOmegaExt~{i}; Jacobian JSur; Integration I1;}
+				Else
+          // dn u_i = i*omega*phi_k on Gamma_i^k
+          For kk In {0:#ListGammaExt~{i}()-1}
+            k=ListGammaExt~{i}(kk);
+            Galerkin{[-I[]*kap*Dof{phi~{i}~{k}}, {u~{i}}];
+              In GammaExt~{i}~{k}; Jacobian JSur; Integration I1;}
+          EndFor
+          // Operator T (Mandatory : Order = 2 + Corner Correction)
+          For kk In {0:#ListGammaExt~{i}()-1}
+            k=ListGammaExt~{i}(kk);
+            Galerkin{[Dof{phi~{i}~{k}}, {phi~{i}~{k}}];
+              In GammaExt~{i}~{k}; Jacobian JSur; Integration I1;}
+            Galerkin{[-Dof{u~{i}}, {phi~{i}~{k}}];
+              In GammaExt~{i}~{k}; Jacobian JSur; Integration I1;}
+            Galerkin{[1/(2*kap^2)*Dof{d phi~{i}~{k}}, {d phi~{i}~{k}}];
+              In GammaExt~{i}~{k}; Jacobian JSur; Integration I1;}
+            If(CORNER_CONDITION == 2)
+              For index In {0:#ListdGammaExt~{i}~{k}()-1:2}
+                j=ListdGammaExt~{i}~{k}(index);
+                l=ListdGammaExt~{i}~{k}(index+1);
+                If (k!=l) // corner point
+                  Galerkin{[-Conj[coef_gamma[]]/2/kap/kap/coef_sigma[]*Dof{phi~{i}~{k}}, {phi~{i}~{k}}];
+                    In Region[{Apt~{i}~{j}~{k}~{l}}]; Jacobian JLin; Integration I1;}
+                  Galerkin{[$TCSource?-h~{i}~{j}~{k}~{l}[]/2/kap/kap:0., {phi~{i}~{k}}];
+                    In Region[{Apt~{i}~{j}~{k}~{l}}]; Jacobian JLin; Integration I1;}
+                //TODO: Flatpoint 
+                //Else 
+                EndIf
+              EndFor
+            EndIf
+          EndFor
+        EndIf
+      }
+    }
+  EndFor
+
+  //Resolution h
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    For kk In {0:#ListGammaExt~{i}()-1}
+      k=ListGammaExt~{i}(kk);
+      For index In {0:#ListdGammaExt~{i}~{k}()-1:2}
+        j=ListdGammaExt~{i}~{k}(index);
+        l=ListdGammaExt~{i}~{k}(index+1);
+        { Name H~{j}~{i}~{l}~{k}; Type FemEquation;
+        	Quantity{
+        	  { Name phi~{i}~{k} ; Type Local; NameOfSpace H_phi~{i}~{k};}
+        	  { Name h~{j}~{i}~{l}~{k} ; Type Local; NameOfSpace H_h~{j}~{i}~{l}~{k};}
+    	    }
+      	  Equation{
+      	    Galerkin{[Dof{h~{j}~{i}~{l}~{k}},{h~{j}~{i}~{l}~{k}}];
+      	      In Apt~{i}~{j}~{k}~{l}; Jacobian JLin; Integration I1;}
+    	      Galerkin{[$TCSource?Conj[coef_sigma[]]/coef_sigma[]*h~{i}~{j}~{k}~{l}[]:0.,{h~{j}~{i}~{l}~{k}}];
+    	        In Apt~{i}~{j}~{k}~{l}; Jacobian JLin; Integration I1;}
+    	      Galerkin{[(Conj[coef_sigma[]*coef_gamma[]]/coef_sigma[]/coef_sigma[] + coef_gamma[]/coef_sigma[]) *{phi~{i}~{k}}, {h~{j}~{i}~{l}~{k}}];
+    	        In Apt~{i}~{j}~{k}~{l}; Jacobian JLin; Integration I1;}
+          }
+        }
+      EndFor
+    EndFor
+  EndFor
+}
+
+Resolution{
+  {Name DDM_2;
+    System{
+      For ii In {0:#myD()-1}
+      	i = myD(ii);
+	      {Name DDM_2~{i}; NameOfFormulation DDM_2~{i}; Type Complex; }
+        For jj In {0:#myD~{i}()-1}
+          j = myD~{i}(jj);
+          {Name G~{j}~{i}; NameOfFormulation G~{j}~{i}; Type Complex; }
+        EndFor
+        If(ABC_ORDER != 0)
+          For kk In {0:#ListGammaExt~{i}()-1}
+            k=ListGammaExt~{i}(kk);
+            For index In {0:#ListdGammaExt~{i}~{k}()-1:2}
+              j=ListdGammaExt~{i}~{k}(index);
+              l=ListdGammaExt~{i}~{k}(index+1);
+              {Name H~{j}~{i}~{l}~{k}; NameOfFormulation H~{j}~{i}~{l}~{k}; Type Complex; }
+            EndFor
+          EndFor
+        EndIf
+      EndFor
+    }
+    Operation{
+      // Initialisation
+      Evaluate[$PhysicalSource = 1];
+      Evaluate[$TCSource = 0];
+      For ii In {0:#myD()-1}
+        i = myD(ii);
+        UpdateConstraint[DDM_2~{i}, dOmegaInt~{i}, Assign];
+        Generate[DDM_2~{i}];
+        Solve[DDM_2~{i}];
+        For jj In {0:#myD~{i}()-1}
+          j = myD~{i}(jj);
+          Generate[G~{j}~{i}];
+          Solve[G~{j}~{i}];
+          PostOperation[G~{j}~{i}];
+        EndFor
+        If(ABC_ORDER != 0)
+          For kk In {0:#ListGammaExt~{i}()-1}
+            k=ListGammaExt~{i}(kk);
+            For index In {0:#ListdGammaExt~{i}~{k}()-1:2}
+              j=ListdGammaExt~{i}~{k}(index);
+              l=ListdGammaExt~{i}~{k}(index+1);    
+              Generate[H~{j}~{i}~{l}~{k}];
+              Solve[H~{j}~{i}~{l}~{k}];
+              PostOperation[H~{j}~{i}~{l}~{k}];
+            EndFor
+          EndFor
+        EndIf
+      EndFor
+      // Set uinc=0 and unlock quantity on Sigma
+      Evaluate[$PhysicalSource = 0];
+      Evaluate[$TCSource = 1];
+      For ii In {0:#myD()-1}
+      	i = myD(ii);
+      	UpdateConstraint[DDM_2~{i}, dOmegaInt~{i}, Assign];
+      EndFor
+
+      IterativeLinearSolver["I-A",SOLVER, TOL, MAXIT, RESTART, ListOfFields_ddm2(), {}, {}]
+      {
+        For ii In {0:#myD()-1}
+          i = myD(ii);
+          GenerateRHSGroup[DDM_2~{i}, Region[{Sigma~{i}, TrdOmegaInt~{i}, dOmegaExt~{i}, dGammaExt}]];
+          SolveAgain[DDM_2~{i}];
+        EndFor
+        For ii In {0:#myD()-1}
+          i = myD(ii);
+          For jj In {0:#myD~{i}()-1}
+            j = myD~{i}(jj);
+            GenerateRHSGroup[G~{j}~{i},Region[{Sigma~{i}~{j}}]];
+            SolveAgain[G~{j}~{i}];
+          EndFor
+          If(ABC_ORDER != 0)
+            For kk In {0:#ListGammaExt~{i}()-1}
+              k=ListGammaExt~{i}(kk);
+              For index In {0:#ListdGammaExt~{i}~{k}()-1:2}
+                j=ListdGammaExt~{i}~{k}(index);
+                l=ListdGammaExt~{i}~{k}(index+1);
+                GenerateRHSGroup[H~{j}~{i}~{l}~{k}, Region[{Apt~{j}~{i}~{l}~{k}}]];
+                SolveAgain[H~{j}~{i}~{l}~{k}];
+              EndFor
+            EndFor
+          EndIf
+        EndFor
+      
+        For ii In {0:#myD()-1}
+          i = myD(ii);
+          For jj In {0:#myD~{i}()-1}
+            j = myD~{i}(jj);
+            PostOperation[G~{j}~{i}];
+          EndFor
+          If(ABC_ORDER != 0)
+            For kk In {0:#ListGammaExt~{i}()-1}
+              k=ListGammaExt~{i}(kk);
+              For index In {0:#ListdGammaExt~{i}~{k}()-1:2}
+                j=ListdGammaExt~{i}~{k}(index);
+                l=ListdGammaExt~{i}~{k}(index+1);
+                PostOperation[H~{j}~{i}~{l}~{k}];
+              EndFor
+            EndFor
+          EndIf
+        EndFor
+      }
+      // Compute solution
+      Evaluate[$PhysicalSource = 1];
+      Evaluate[$TCSource = 1];
+
+      For ii In {0:#myD()-1}
+        i = myD(ii);
+        UpdateConstraint[DDM_2~{i}, dOmegaInt~{i}, Assign];
+        GenerateRHSGroup[DDM_2~{i}, Region[{Sigma~{i}, TrdOmegaInt~{i}, dOmegaExt~{i}, dGammaExt}]];
+        SolveAgain[DDM_2~{i}];
+      EndFor
+      PostOperation[DDM];
+    }
+  }
+}
+
+// Post Processing and Post Operation specific to DDM-2 algorithm
+// ================
+PostProcessing{
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    For kk In {0:#ListGammaExt~{i}()-1}
+      k=ListGammaExt~{i}(kk);
+      For index In {0:#ListdGammaExt~{i}~{k}()-1:2}
+        j=ListdGammaExt~{i}~{k}(index);
+        l=ListdGammaExt~{i}~{k}(index+1);
+        {Name H~{j}~{i}~{l}~{k}; NameOfFormulation H~{j}~{i}~{l}~{k};
+          Quantity {
+           {Name h~{j}~{i}~{l}~{k}; Value {Local { [{h~{j}~{i}~{l}~{k}}] ; In Apt~{j}~{i}~{l}~{k}; Jacobian JLin; }}}
+          }
+        }
+      EndFor
+    EndFor
+  EndFor
+}
+
+PostOperation{
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    For kk In {0:#ListGammaExt~{i}()-1}
+      k=ListGammaExt~{i}(kk);
+      For index In {0:#ListdGammaExt~{i}~{k}()-1:2}
+        j=ListdGammaExt~{i}~{k}(index);
+        l=ListdGammaExt~{i}~{k}(index+1);
+        {Name H~{j}~{i}~{l}~{k}; NameOfPostProcessing H~{j}~{i}~{l}~{k};
+          Operation {Print [h~{j}~{i}~{l}~{k}, OnElementsOf Apt~{j}~{i}~{l}~{k}, StoreInField tag_h~{j}~{i}~{l}~{k}];}
+        }
+      EndFor
+    EndFor
+  EndFor
+}

DDM-Corner-Helmholtz-2D/src/ABC_Decomposition.pro

+// Dispatch subdomain to (MPI) processes (list myD())
+// Provide a unique tag to every transmitted data (e.g.: g_{i,j})
+Function{
+  myD = {} ; // the domains current process will treat (if MPI_NUM_PROC = 1 then myD()= {0,1,2,..., NDOM-1}
+  For i In {0:(NDom-1)}
+    myD~{i} = {};
+    If (i % MPI_Size == MPI_Rank)
+      myD() += D(i);
+      myD~{i} += D~{i}();
+    EndIf
+  EndFor
+  // Field managed by current MPI Process
+  // If sequential then every field
+  ListOfFields_g = {};
+  ListOfConnectedFields_g = {};
+  ListOfFields_h = {};
+  ListOfConnectedFields_h = {};
+
+  // Loop on subdomains
+  For ii In {0:#myD()-1}
+    i = myD(ii);
+    //Neighbors
+    For jj In {0:#myD~{i}()-1}
+      j = myD~{i}(jj);
+
+      // Quantities coming from "the other side"
+      tag_g~{i}~{j} = 10000 + i*50 + j;
+      ListOfFields_g() += tag_g~{i}~{j};
+      g~{i}~{j}[ Sigma~{i}~{j} ] = ComplexScalarField[XYZ[]]{ tag_g~{i}~{j} };
+      // Quantities to be transmitted to domain j
+      tag_g~{j}~{i} = 10000 + j*50 + i;
+      ListOfConnectedFields_g() += 1;
+      ListOfConnectedFields_g() += tag_g~{j}~{i};
+
+      // RHS (for jacobi)
+      tag_b~{i}~{j} = 11000 + i*50 + j;
+      tag_b~{j}~{i} = 11000 + j*50 + i;
+      b~{i}~{j}[ Sigma~{i}~{j} ] = ComplexScalarField[XYZ[]]{ tag_b~{i}~{j} };
+    EndFor
+    // For DDM-2
+    For kk In {0:#ListGammaExt~{i}()-1}
+      k=ListGammaExt~{i}(kk);
+      For index In {0:#ListdGammaExt~{i}~{k}()-1:2}
+        j=ListdGammaExt~{i}~{k}(index);
+        l=ListdGammaExt~{i}~{k}(index+1);
+        // Quantities coming from "the other side" of Aijkl
+        tag_h~{i}~{j}~{k}~{l} = 100000 + i*1000 + j*100 + 10*k + l;
+        ListOfFields_h() += tag_h~{i}~{j}~{k}~{l};
+        h~{i}~{j}~{k}~{l}[ Apt~{i}~{j}~{k}~{l} ] = ComplexScalarField[XYZ[]]{ tag_h~{i}~{j}~{k}~{l} };
+        // Quantities to be transmitted to domain j
+        tag_h~{j}~{i}~{l}~{k} = 100000 + j*1000 + i*100 + 10*l + k;
+        ListOfConnectedFields_h() += 1;
+        ListOfConnectedFields_h() += tag_h~{j}~{i}~{l}~{k};
+      EndFor
+    EndFor
+  EndFor
+}

DDM-Corner-Helmholtz-2D/src/ABC_MonoDomain.pro

+// Mono domain problem
+
+Constraint {
+  {Name Dirichlet; Type Assign; Case{{Region dOmegaInt; Value -uInc[]; } } }
+}
+
+FunctionSpace{
+  {Name H_grad; Type Form0;
+    BasisFunction{
+       {Name u; NameOfCoef ui; Function BF_Node;
+	      Support Region[{Omega, dOmegaExt, dOmegaInt}]; Entity NodesOf[All];}
+      }
+    Constraint{
+      {NameOfCoef ui; EntityType NodesOf; NameOfConstraint Dirichlet;}
+    }
+  }
+  For k In {0:(#ListGammaExt()-1)}
+    {Name H_phi~{k}; Type Form0;
+      BasisFunction{
+	      {Name phi~{k}; NameOfCoef phik; Function BF_Node;
+	      Support Region[{GammaExt~{k}, dGammaExt~{k}}]; Entity NodesOf[All];}
+      }
+    }
+  EndFor
+}
+
+Formulation {
+  // Formulation for a Dirichlet boundary condition
+  { Name MonoDomain; Type FemEquation;
+    Quantity{
+      { Name u ; Type Local; NameOfSpace H_grad;}
+        For k In {0:(#ListGammaExt()-1)}
+          { Name phi~{k} ; Type Local; NameOfSpace H_phi~{k};}
+        EndFor
+      }
+    Equation{
+      //MonoDomain equation
+      Galerkin{[Dof{Grad u}, {Grad u}];
+        In Omega; Jacobian JVol; Integration I1;}
+      Galerkin{[-kap^2*Dof{u}, {u}];
+        In Omega; Jacobian JVol; Integration I1;}
+      // ABC
+      If(ABC_ORDER == 0)
+        Galerkin{[-I[]*kap*Dof{u}, {u}];
+          In dOmegaExt; Jacobian JSur; Integration I1;}
+      ElseIf(ABC_ORDER == 2)
+        For k In {0:(#ListGammaExt()-1)}
+          Galerkin{[-I[]*kap*Dof{phi~{k}}, {u}];
+            In GammaExt~{k}; Jacobian JSur; Integration I1;}
+          Galerkin{[Dof{phi~{k}}, {phi~{k}}];
+            In GammaExt~{k}; Jacobian JSur; Integration I1;}
+          Galerkin{[-Dof{u}, {phi~{k}}];
+            In GammaExt~{k}; Jacobian JSur; Integration I1;}
+          Galerkin{[1/(2*kap^2)*Dof{d phi~{k}}, {d phi~{k}}];
+              In GammaExt~{k}; Jacobian JSur; Integration I1;}
+          If(CORNER_CONDITION == 2)
+            For index In {0:#ListdGammaExt~{k}()-1}
+              // Loop on end-point
+							l=ListdGammaExt~{k}(index); // Index of ajacent segment
+              Galerkin{[-I[]/(4*kap)*(Cos[theta]/Cos[theta/2]+Cos[theta/2])*Dof{phi~{k}}, {phi~{k}}];
+                In Apt~{k}~{l}; Jacobian JLin; Integration I1;}
+              Galerkin{[-I[]/(4*kap)*(-Cos[theta]/Cos[theta/2]+Cos[theta/2])*Dof{phi~{l}}, {phi~{k}}];
+                In Apt~{k}~{l}; Jacobian JLin; Integration I1;}
+            EndFor
+          EndIf
+        EndFor
+      EndIf
+    }
+  }
+}
+
+Resolution{
+  {Name MonoDomain;
+    System{ {Name A; NameOfFormulation MonoDomain; Type Complex; } }
+    Operation{
+      CreateDir[ABC_DIR];
+      Generate[A];
+      Solve[A];
+      PostOperation[MonoDomain];
+    }
+  }
+}
+
+PostProcessing{
+  {Name MonoDomain; NameOfFormulation MonoDomain;
+    Quantity {
+      {Name u; Value {Local { [{u}] ; In Omega; Jacobian JVol; }}}
+      {Name uNorm; Value {Local { [Norm[{u}]] ; In Omega; Jacobian JVol; }}}
+      {Name err; Value {Local { [Norm[{u}-uex[]]] ; In Omega; Jacobian JVol; }}}
+    }
+  }
+}
+
+PostOperation{
+  {Name MonoDomain; NameOfPostProcessing MonoDomain;
+    Operation {
+      Print [u, OnElementsOf Omega, Name "Scattered_Field", File StrCat[ABC_DIR, StrCat[ABC_PREFIX, "u_mono.pos"]]];
+      Print [uNorm, OnElementsOf Omega, Name "Scattered_Field_(abs._value)", File StrCat[ABC_DIR,StrCat[ABC_PREFIX, "u_mono_abs.pos"]]];
+      Print [err, OnElementsOf Omega, Name "Error_numeric_vs_analytic", File StrCat[ABC_DIR,StrCat[ABC_PREFIX, "u_mono_err.pos"]]];
+    }
+  }
+}

DDM-Corner-Helmholtz-2D/src/ABC_main.data

+DefineConstant[
+	//Geometry
+	theta = {-2*Pi/3, ReadOnly 1, Name "Input/Geometry/0Theta"},
+	NGammaExt = {6, Visible 0, Name "Input/Geometry/1Nb. segments"},
+	r = {5, Min 0, Max 100, Step 0.1, Name "Input/Geometry/2Radius of the ABC"},
+	rs = {1, Min 1, Max 10, Step 0.1, Name "Input/Geometry/4Radius of the obstacle"},
+	xs = {0, Min -1, Max 1, Step 0.1, Name "Input/Geometry/6X-location of the obstacle"},
+	ys = {0, Min -1, Max 1, Step 0.1, Name "Input/Geometry/8Y-location of the obstacle"},
+	//Input
+	kap = {2*Pi, Min 1, Max 30, Step 0.1, Name "Input/Input/0Wavenumber"},
+	// Mesh
+	NLambda = {15, Min 10, Max 100, Step 1, Name "Input/Mesh/NLambda"},
+	h = {2*Pi/NLambda/kap, ReadOnly 1, Name "Input/Mesh/h"},
+	// DDM
+	NDom = {3, Choices {2="2", 3="3", 6="6"}, Name "Input/DDM/6Nb. sub-domains"}
+];
+

DDM-Corner-Helmholtz-2D/src/ABC_main.geo

+//Hexagone regulier et point source/obstacle disque
+Mesh.CharacteristicLengthMin = h;
+Mesh.CharacteristicLengthMax = h;
+
+Point(10) = {xs,ys,0,h}; //Center of the disk and hexagone
+For i In {0:5}
+	Point(i+1) = {r*Cos[i*Pi/3],r*Sin[i*Pi/3],0,h}; //6 points hexagone
+	Point(10+i+1) = {xs+rs*Cos[i*Pi/3],ys+rs*Sin[i*Pi/3],0,h}; //6 points disk (obstacle)
+EndFor
+For i In {1:6}
+	Line(i) = {i,(i%6)+1}; //segments of the hexagone
+	Circle(10+i) = {10+i,10,10+(i%6)+1}; //circle arc of the disk
+EndFor
+
+// Segments of the hexagone (numbering does not depend on the decomposition)
+For k In {0:(NGammaExt-1)}
+	Physical Curve(30+k) = {k+1}; //segments Gamma k ext
+	Physical Point(40+k) = {k+1}; //noeud k-1,k
+	Physical Point(50+k) = {10+k+1}; //Noeuds Gamma Int
+EndFor
+
+// Numbering depends on the number of subdomains
+If (NDom==1)
+	Curve Loop(1) = {1:6}; //hexagone
+	Curve Loop(2) = {11:16}; // disk
+	Plane Surface(1) = {1,2}; // Hexagone minus disk (=propagation domain)
+
+	Physical Surface(0) = {1}; //Propagation domaine
+	Physical Curve(10) = {11:16}; // Boundary of disk (Gamma Int)
+EndIf
+
+// With 2 subdomaines, Omega_0 and Omega_1 share two boundaries in common
+If (NDom==2)
+	Line(21) = {1,11};
+	Line(22) = {4,14};
+	Curve Loop(1) = {1,2,3,22,-13,-12,-11,-21}; //Boundary of Omega 0
+	Plane Surface(1) = {1}; // Omega 0
+	Curve Loop(2) = {21,-16,-15,-14,-22,4,5,6}; //Boundary of Omega 1
+	Plane Surface(2) = {2}; //Omega 1
+
+	Physical Surface(0) = {1}; //Omega 0
+	Physical Surface(1) = {2}; //Omega 1
+
+	Physical Curve(10) = {11,12,13}; //Gamma Int 0
+	Physical Curve(11) = {14,15,16}; //Gamma Int 1
+
+	Physical Curve(20) = {21,22}; //Sigma_{0,1}
+EndIf
+
+If (NDom==3)
+	For i In {0:(NDom-1)}
+		Line(21+i) = {1+i*(NGammaExt/NDom),11+i*(NGammaExt/NDom)};
+	EndFor
+
+	Curve Loop(1) = {1,2,22,-12,-11,-21}; // Boundary of Omega 0
+	Plane Surface(1) = {1}; //Omega 0
+
+	Curve Loop(2) = {3,4,23,-14,-13,-22};
+	Plane Surface(2) = {2};
+
+	Curve Loop(3) = {21,-16,-15,-23,5,6};
+	Plane Surface(3) = {3};
+
+	For i In {0:(NDom-1)}
+		Physical Surface(i) = {i+1}; //Omega_i
+		Physical Curve(20+i) = {20+i+1}; //Sigma_{i, (i-1)}
+		Physical Curve(10+i) = {11 + 2*i, 12 + 2*i}; //Gamma Int of Omega_i
+	EndFor
+EndIf
+
+If (NDom==6)
+	For j In {1:NDom}
+		Line(20+j) = {j,10+j}; // Sigma
+	EndFor
+
+	For i In {1:(NDom-1)}
+		Curve Loop(i) = {i,20+i+1,-10-i,-20-i}; //Boundary Omega_i
+		Plane Surface(i) = {i}; //Omega_i
+	EndFor
+	Curve Loop(6) = {6,21,-16,-26}; //Boundary Omega_5
+	Plane Surface(6) = {6}; //Omega_5
+
+	For i In {0:(NDom-1)}
+		Physical Surface(i) = {i+1}; //Omega_i
+		Physical Curve(10+i) = {10+i+1}; //Gamma Interior of Omega_i
+		Physical Curve(20+i) = {20+i+1}; //Sigma_{i-1, i}		
+	EndFor
+EndIf
+
+Solver.AutoMesh = 2; // auto mesh

DDM-Corner-Helmholtz-2D/src/ABC_main.pro

+DefineConstant[
+	//Input
+	alpha = {0, Min 0, Max 2*Pi, Step 0.1, Name "Input/Input/5Incident angle"},
+	// ABC
+	ABC_ORDER = {2, Choices {0="Order 0", 2="Order 2"}, Name "Input/ABC/0ABC order condition"},
+	CORNER_CONDITION = {2, Choices {1 ="Homog. Neumann", 2="Corner Correction"}, Name "Input/ABC/3Corner Condition", ReadOnly (ABC_ORDER == 0)},
+	// DDM
+	beta = {-kap*Cos[theta], ReadOnly 1, Name "Input/DDM/8beta parameter"},
+	// Solver
+	SOLVER = {"jacobi", Choices {"jacobi", "gmres", "print"}, Name "Input/IterativeSolver/0Solver"},
+	TOLlog10 = {-6, Max -1, Min -16, Step -1, Name "Input/IterativeSolver/1Tolerance (log10)"},
+	TOL = 10^(TOLlog10),
+	MAXIT = {10, Min 10, Max 10000, Step 1, Name "Input/IterativeSolver/2Max it"},
+	RESTART = {0, Min 0, Max 1000, Step 1, Name "Input/IterativeSolver/3Restart"},
+		//Output
+	ABC_DIRREL = {"out_ABC/", Name "Input/Output/01Output Directory"},
+	ABC_DIR = {StrCat[WorkingDir, ABC_DIRREL], ReadOnly 1, Visible 0, Name "Input/Output/Output Directory (Absolute)"}	
+	ABC_PREFIX = {"", Name "Input/Output/05Prefix for filename"}
+];
+
+Function {
+  I[] = Complex[0.,1.]; // Pure imaginary "i"
+  alpha_vect[] = Vector[Cos[alpha],Sin[alpha],0]; // incidence vector
+  uInc[] = Complex[Cos[kap*(XYZ[]*alpha_vect[])],Sin[kap*(XYZ[]*alpha_vect[])]]; // Incident plane wave
+  //Exact solution for a cylindrical obstacle located at (0,0,0) of radius = rs and incident angle = alpha
+  uex[] = AcousticFieldSoftCylinder[XYZ[]]{kap, rs, alpha};
+}
+
+// Mono domain case
+Group{
+  Omega = Region[{(0:(NDom-1))}];
+  GammaExt = Region[{30:35}];
+  Sigma = Region[{20:25}];
+  dOmegaInt = Region[{10:15}];
+  dOmegaExt = Region[{GammaExt}];
+  ListGammaExt() = {}; // indices k of the segments
+  // Fourier or Neumann condition
+  For k In {0:NGammaExt-1}
+    ListGammaExt() +=k;
+    GammaExt~{k} = Region[(30+k)];
+    // End-points
+    l_left = (k-1+NGammaExt)%NGammaExt;
+    l_right = (k+1)%NGammaExt;
+    ListdGammaExt~{k}() = {l_left,l_right};
+    Apt~{k}~{l_left} = Region[(40 + k)];
+    Apt~{k}~{l_right} = Region[(40 + (k+1)%6)];
+    // Every end-points of Gamma_k
+    dGammaExt~{k} = Region[{Apt~{k}~{l_left},Apt~{k}~{l_right}}];
+  EndFor
+  dGammaExt = Region[{40:45}]; //every end-points
+  TrdOmegaInt = ElementsOf[Omega, OnOneSideOf dOmegaInt];
+}
+
+//DDM Case
+If (NDom == 2)
+  Group {
+  D() = {};
+  For i In {0:(NDom-1)}
+    Omega~{i} = Region[{i}];
+    D() += i;
+    // left = right = {1, 0}
+    j = (i+1)%NDom; // The other subdomain
+    D~{i} = {j};
+
+    GammaInt~{i} = Region[{(10+i)}];
+    NGammaExt~{i} = (NGammaExt/NDom); // Number of segments Gamma^i_k = dOmega_i cap Gamma
+    ListGammaExt~{i}() = {}; //Set of indices k such that Gamma^i_k is not empty
+    // Every subdomain has 3 Gamma^i_k :  
+    // if i == 0 then k = {0,1,2} 
+    // else i == 1 then k = {3,4,5} 
+    For kk In {0:(NGammaExt~{i}-1)} // kk = 0,1,2
+      k = kk + 3*i; // k = 0,1,2 or 3,4,5
+      ListGammaExt~{i}() += k;
+      ListdGammaExt~{i}~{k}() = {}; // set of indices (j,l) of every end-points of Gamma^i_k
+      GammaExt~{i}~{k} = Region[{(30+k)}]; //Gamma^i_k
+      // Subdomain number of neightbor segment
+      If(kk == 0)
+        j_left = j;
+      Else
+        j_left = i;
+      EndIf
+      If(kk == 2)
+        j_right = j;
+      Else
+        j_right = i;
+      EndIf
+      // Indice of the previous/next segment (independant of subdomain number)
+      l_left = (k-1+NGammaExt)%NGammaExt;
+      l_right = (k+1)%NGammaExt;
+ 
+      // first end-point A^ij_kl of Gamma^i_k
+      Apt~{i}~{j_left}~{k}~{l_left} = Region[{(40+k)}];
+      dGammaExt~{i}~{j_left}~{k}~{l_left} = Region[{(40+k)}];
+      ListdGammaExt~{i}~{k}() += j_left; // j
+      ListdGammaExt~{i}~{k}() += l_left; // l
+
+      //second end-point A^ij_kl of Gamma^i_k
+      Apt~{i}~{j_right}~{k}~{l_right} = Region[{(40+(k+1)%NGammaExt)}];
+      dGammaExt~{i}~{j_right}~{k}~{l_right} = Region[{(40+(k+1)%NGammaExt)}];
+      ListdGammaExt~{i}~{k}() += j_right; // j
+      ListdGammaExt~{i}~{k}() += l_right; // l
+
+      //Set C^i_k of every end-points of Gamma^i_k
+      Cset~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
+      dGammaExt~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
+    EndFor
+
+    dOmegaInt~{i}= Region[{GammaInt~{i}}];
+    dOmegaExt~{i}= Region[{GammaExt~{i}~{3*i}, GammaExt~{i}~{1 + 3*i}, GammaExt~{i}~{2 + 3*i}}];
+    //Sigma_{i,j} = dOmega_i cap dOmega_j
+    Sigma~{i}~{j} = Region[{20}];
+
+    //2 end points on Sigma cap GammaExt
+    ListdSigmaExt~{i}~{j}() = {};
+    //First
+    k = 0 + 3*i;
+    l_left = (k - 1 + NGammaExt)%NGammaExt;
+    dSigmaExt~{i}~{j}~{0} = Region[{Apt~{i}~{j}~{k}~{l_left}}];
+    ListdSigmaExt~{i}~{j}() += j;
+    ListdSigmaExt~{i}~{j}() += l_left;
+    //Second
+    k = 2 + 3*i;
+    l_right = (k +1 )%NGammaExt;
+    dSigmaExt~{i}~{j}~{1} = Region[{Apt~{i}~{j}~{k}~{l_right}}];
+    ListdSigmaExt~{i}~{j}() += j;
+    ListdSigmaExt~{i}~{j}() += l_right;
+    // union
+    dSigmaExt~{i}~{j} = Region[{dSigmaExt~{i}~{j}~{0}, dSigmaExt~{i}~{j}~{1}}];
+
+    //2 end points on GammaInt
+    l = j;
+    ListdSigmaInt~{i}~{j}() = {};
+    dSigmaInt~{i}~{j}~{0} = Region[{(50)}];
+    ListdSigmaInt~{i}~{j}() += i; // k == i (6 sub dom)
+    ListdSigmaInt~{i}~{j}() += l; // l == j (6 sub dom)
+    dSigmaInt~{i}~{j}~{1} = Region[{(53)}];
+    ListdSigmaInt~{i}~{j}() += i; // k == i (6 sub dom)
+    ListdSigmaInt~{i}~{j}() += l; // l == j (6 sub dom)
+    dSigmaInt~{i}~{j} = Region[{dSigmaInt~{i}~{j}~{0}, dSigmaInt~{i}~{j}~{1}}];
+
+    dSigma~{i}~{j} = Region[{dSigmaExt~{i}~{j},dSigmaInt~{i}~{j}}];
+
+    Sigma~{i} = Region[{Sigma~{i}~{j}}];
+  EndFor
+  }
+EndIf
+
+If (NDom==3)
+  Group {
+  D() = {};
+  For i In {0:(NDom-1)}
+    Omega~{i} = Region[{i}];
+    D() += i;
+    left = (i-1+NDom)%NDom; //previous subdomain (counter-clowise orientation)
+    right = (i+1)%NDom;  //next subdomain (counter-clowise orientation)
+    D~{i} = {left, right}; // neighbors subdomains
+
+    dOmegaInt~{i} = Region[{(10+i)}];
+    NGammaExt~{i} = (NGammaExt/NDom); // Number of Gamma^i_k (= 2)
+    ListGammaExt~{i}() = {}; //Indices k of these Gamma^i_k
+
+    For kk In {0:(NGammaExt~{i}-1)} // kk = 0,1
+      k = kk + 2*i; // k = [0,1] or [2,3] or [4,5]
+      ListGammaExt~{i}() += k;
+      ListdGammaExt~{i}~{k}() = {}; // indices (j,l) of the neightbor Gamma^i_k
+      GammaExt~{i}~{k} = Region[{(30+k)}]; //Gamma^i_k
+      // Subdomain number of neightbor segment
+      If(kk == 0)
+        j_left = left;
+        j_right = i;
+      EndIf
+      If(kk == 1)
+        j_left = i;
+        j_right = right;
+      EndIf
+      l_left = (k-1+NGammaExt)%NGammaExt;
+      l_right = (k+1)%NGammaExt;
+
+      //first end-point A^ij_kl of Gamma^i_k
+      Apt~{i}~{j_left}~{k}~{l_left} = Region[{(40+k)}];
+      ListdGammaExt~{i}~{k}() += j_left; // j
+      ListdGammaExt~{i}~{k}() += l_left; // l
+
+      //second noeud A^ij_kl de Gamma^i_k
+      Apt~{i}~{j_right}~{k}~{l_right} = Region[{(40+(k+1)%NGammaExt)}];
+      ListdGammaExt~{i}~{k}() += j_right; // j
+      ListdGammaExt~{i}~{k}() += l_right; // l
+
+      //ensemble C^i_k = noeuds bord Gamma^i_k
+      Cset~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
+      dGammaExt~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
+    EndFor
+    // Union
+    dOmegaExt~{i} = Region[{GammaExt~{i}~{2*i}, GammaExt~{i}~{1 + 2*i}}];
+
+    // Transmission boundaries
+    Sigma~{i}~{left} = Region[{(20 + i)}];
+    Sigma~{i}~{right} = Region[{(20 + right)}];
+    Sigma~{i} = Region[{Sigma~{i}~{left}, Sigma~{i}~{right}}];
+    // Unique End-point of Sigma_{ij} that also belons to Gamma
+    ListdSigmaExt~{i}~{left}() = {};
+    ListdSigmaExt~{i}~{right}() = {};
+    k_left = 2*i;
+    l_left = (k_left - 1 +NGammaExt )%NGammaExt;
+    k_right = 2*i +1 ;
+    l_right = (k_right + 1 )%NGammaExt;
+    dSigmaExt~{i}~{left}~{0} = Region[{Apt~{i}~{left}~{k_left}~{l_left}}];
+    ListdSigmaExt~{i}~{left}() += k_left;
+    ListdSigmaExt~{i}~{left}() += l_left;
+    dSigmaExt~{i}~{left} = Region[{dSigmaExt~{i}~{left}~{0}}];
+
+    dSigmaExt~{i}~{right}~{0} = Region[{Apt~{i}~{right}~{k_right}~{l_right}}];
+    ListdSigmaExt~{i}~{right}() += k_right;
+    ListdSigmaExt~{i}~{right}() += l_right;
+    dSigmaExt~{i}~{right} = Region[{dSigmaExt~{i}~{right}~{0}}];
+
+    ListdSigmaInt~{i}~{left}() = {};
+    dSigmaInt~{i}~{left}~{0} = Region[{(50 + 2*i)}];
+    ListdSigmaInt~{i}~{left}() += i; // Useless ?
+    ListdSigmaInt~{i}~{left}() += left; // Useless ?
+    dSigmaInt~{i}~{left} = Region[{dSigmaInt~{i}~{left}~{0}}];
+
+    ListdSigmaInt~{i}~{right}() = {};
+    dSigmaInt~{i}~{right}~{0} = Region[{(50 + 2*i + 1)}];
+    ListdSigmaInt~{i}~{right}() += i; // Useless ?
+    ListdSigmaInt~{i}~{right}() += right; // Useless ?
+    dSigmaInt~{i}~{right} = Region[{dSigmaInt~{i}~{right}~{0}}];
+
+    dSigma~{i}~{left} = Region[{dSigmaExt~{i}~{left},dSigmaInt~{i}~{left}}];
+    dSigma~{i}~{right} = Region[{dSigmaExt~{i}~{right},dSigmaInt~{i}~{right}}];
+
+    //Sigma i = Sigma i i-1 U Sigma i i+1
+    Sigma~{i} = Region[{Sigma~{i}~{left},Sigma~{i}~{right}}];
+    dSigmaExt~{i} = Region[{dSigmaExt~{i}~{left},dSigmaExt~{i}~{right}}];
+    dSigmaInt~{i} = Region[{dSigmaInt~{i}~{left},dSigmaInt~{i}~{right}}];
+    dSigma~{i} = Region[{dSigmaExt~{i},dSigmaInt~{i}}];
+  EndFor
+  }
+EndIf
+
+If (NDom==6)
+  Group {
+  D() = {};
+  For i In {0:(NDom-1)}
+    Omega~{i} = Region[{i}]; //Omega_i
+    D() += i;
+    left = (i-1+NDom)%NDom; // previous subdomain (counter-clockwise orientation)
+    right = (i+1)%NDom; //next subdomain (counter-clockwise orientation)
+    D~{i} = {left,right};
+
+    NGammaExt~{i} = (NGammaExt/NDom); // Number of segments Gamma^i_k (=1)
+    dOmegaInt~{i} = Region[{(10+i)}]; //Gamma i int
+    ListGammaExt~{i}() = {}; // index k of this Gamma^i_k and in this case k=i
+    k = i;
+    ListGammaExt~{i}() += k;
+    GammaExt~{i}~{k} = Region[{(30+k)}]; //Gamma i k ext
+    ListdGammaExt~{i}~{k}() = {}; // list of indices (j,l) for end-points A^{ij}_{kl} of Gamma^i_k
+    j_left = left;
+    l_left = left; // l = j
+    j_right = right; 
+    l_right = right; //l =j
+    // first end-point A^ij_kl of Gamma^i_k
+    Apt~{i}~{j_left}~{k}~{l_left} = Region[{(40+k)}];
+    dGammaExt~{i}~{j_left}~{k}~{l_left} = Region[{(40+k)}];
+    ListdGammaExt~{i}~{k}() += j_left; // j
+    ListdGammaExt~{i}~{k}() += l_left; // l
+    //second end-point A^ij_kl of Gamma^i_k
+    Apt~{i}~{j_right}~{k}~{l_right} = Region[{(40+(k+1)%NGammaExt)}];
+    dGammaExt~{i}~{j_right}~{k}~{l_right} = Region[{(40+(k+1)%NGammaExt)}];
+    ListdGammaExt~{i}~{k}() += j_right; // j
+    ListdGammaExt~{i}~{k}() += l_right; // l
+    //Set C^i_k = end-points of Gamma^i_k
+    Cset~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
+    dGammaExt~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
+    // Union (here =) of gamma^i_k
+    dOmegaExt~{i}= Region[{GammaExt~{i}~{k}}];
+
+    //Sigma_{ij} = dOmega_i cap dOmega_j
+    // For the left subdomain
+    Sigma~{i}~{left} = Region[{(20+i)}];
+    ListdSigmaExt~{i}~{left}() = {}; // get k and l for provided i and j
+    dSigmaExt~{i}~{left}~{0} = Region[{Apt~{i}~{left}~{i}~{left}}];
+    ListdSigmaExt~{i}~{left}() += k; // k == i (6 sub dom)
+    ListdSigmaExt~{i}~{left}() += l_left; // l == left (6 sub dom)
+    dSigmaExt~{i}~{left} = Region[{dSigmaExt~{i}~{left}~{0}}];
+
+    ListdSigmaInt~{i}~{left}() = {};  // get k and l for provided i and j
+    dSigmaInt~{i}~{left}~{0} = Region[{(50+i)}];
+    ListdSigmaInt~{i}~{left}() += i; // k == i (6 sub dom)
+    ListdSigmaInt~{i}~{left}() += left; // l == left (6 sub dom)
+    dSigmaInt~{i}~{left} = Region[{dSigmaInt~{i}~{left}~{0}}];
+
+    dSigma~{i}~{left} = Region[{dSigmaExt~{i}~{left},dSigmaInt~{i}~{left}}];
+
+    // For the right subdomain
+    Sigma~{i}~{right} = Region[{(20+right)}];
+
+    ListdSigmaExt~{i}~{right}() = {};  // get k and l for provided i and j
+    dSigmaExt~{i}~{right}~{0} = Region[{Apt~{i}~{right}~{i}~{right}}];
+    ListdSigmaExt~{i}~{right}() += i; // k == i (6 sub dom)
+    ListdSigmaExt~{i}~{right}() += right; // l == right (6 sub dom)
+    dSigmaExt~{i}~{right} = Region[{dSigmaExt~{i}~{right}~{0}}];
+
+    ListdSigmaInt~{i}~{right}() = {};  // get k and l for provided i and jerer k et l pour i et j donnes
+    dSigmaInt~{i}~{right}~{0} = Region[{(50+(i+1)%NDom)}];
+    ListdSigmaInt~{i}~{right}() += i; // k == i (6 sub dom)
+    ListdSigmaInt~{i}~{right}() += right; // l == right (6 sub dom)
+    dSigmaInt~{i}~{right} = Region[{dSigmaInt~{i}~{right}~{0}}];
+
+    dSigma~{i}~{right} = Region[{dSigmaExt~{i}~{right},dSigmaInt~{i}~{right}}];
+
+    //Sigma i = Sigma_{i, left} U Sigma_{i,right}
+    Sigma~{i} = Region[{Sigma~{i}~{left},Sigma~{i}~{right}}];
+    dSigmaExt~{i} = Region[{dSigmaExt~{i}~{left},dSigmaExt~{i}~{right}}];
+    dSigmaInt~{i} = Region[{dSigmaInt~{i}~{left},dSigmaInt~{i}~{right}}];
+    dSigma~{i} = Region[{dSigmaExt~{i},dSigmaInt~{i}}];
+  EndFor
+  }
+EndIf
+
+
+//===========================
+Include "ABC_MonoDomain.pro";
+Include "ABC_DDM.pro";
+//===========================

DDM-Corner-Helmholtz-2D/src/Common.pro

+// Common quantities for every resoltion
+Jacobian {
+  { Name JVol ; Case { { Region All ; Jacobian Vol ; } } }
+  { Name JSur ; Case { { Region All ; Jacobian Sur ; } } }
+  { Name JLin ; Case { { Region All ; Jacobian Lin ; } } }
+}
+
+Integration {
+  { Name I1 ;
+    Case {
+      { Type Gauss ;
+      	Case {
+	        { GeoElement Point ; NumberOfPoints  1 ; }
+	        { GeoElement Line ; NumberOfPoints  4 ; }
+	        { GeoElement Triangle ; NumberOfPoints  6 ; }
+	        { GeoElement Quadrangle ; NumberOfPoints 7 ; }
+	        { GeoElement Tetrahedron ; NumberOfPoints 15 ; }
+	        { GeoElement Hexahedron ; NumberOfPoints 34 ; }
+      	}
+      }
+    }
+  }
+}
\ No newline at end of file