From 22fe387843381c08461833c5f07c16caead4be63 Mon Sep 17 00:00:00 2001
From: pmarchne <philippe.marchner@siemens.com>
Date: Thu, 8 Apr 2021 17:11:08 +0200
Subject: [PATCH] draft reorganization of helmholtz examples

---
 contrib/bessel/CMakeLists.txt                 |   2 +-
 .../{ => diskScattering}/CMakeLists.txt       |   0
 .../{ => diskScattering}/README.md            |   0
 .../helmholtz2d/{ => diskScattering}/main.cpp |   0
 examples/helmholtz2d/waveguide/Readme.md      |  99 +++
 .../waveguide/longitudinal/CMakeLists.txt     |   8 +
 .../waveguide/longitudinal}/main.cpp          |   1 -
 .../waveguide/longitudinal/main_length.cpp}   |   1 -
 .../waveguide/longitudinal/main_omega.cpp}    |   1 -
 .../waveguide/transverse/CMakeLists.txt       |   8 +
 .../waveguide/transverse}/main.cpp            |   0
 .../waveguide/transverse/main_length.cpp}     |   0
 .../waveguide/transverse/main_omega.cpp}      |   0
 .../flowPastCylinder/CMakeLists.txt           |   8 +
 .../main.cpp                                  |   0
 .../helmholtzFlow/freefield/CMakeLists.txt    |   8 +
 examples/helmholtzFlow/freefield/main.cpp     | 564 ++++++++++++++++++
 .../runTests_circle.sh                        |   8 +-
 .../runTests_square.sh                        |  10 +-
 .../helmholtzDuctHeteroX/CMakeLists.txt       |   8 -
 .../helmholtzDuctHeteroY/CMakeLists.txt       |   8 -
 .../helmholtzflow/Duct_HeteroX/CMakeLists.txt |   8 -
 .../FlowPastCylinder/CMakeLists.txt           |   8 -
 examples/helmholtzFlow/waveguide/Readme.md    |  78 +++
 .../waveguide/nonUniform/CMakeLists.txt       |   8 +
 .../nonUniform}/main.cpp                      |   3 +-
 .../nonUniform/main_omega.cpp}                |   0
 .../waveguide/uniform/CMakeLists.txt          |   8 +
 .../helmholtzFlow/waveguide/uniform/main.cpp  | 266 +++++++++
 .../uniform}/main_omega.cpp                   |   0
 30 files changed, 1066 insertions(+), 47 deletions(-)
 rename examples/helmholtz2d/{ => diskScattering}/CMakeLists.txt (100%)
 rename examples/helmholtz2d/{ => diskScattering}/README.md (100%)
 rename examples/helmholtz2d/{ => diskScattering}/main.cpp (100%)
 mode change 100755 => 100644
 create mode 100644 examples/helmholtz2d/waveguide/Readme.md
 create mode 100644 examples/helmholtz2d/waveguide/longitudinal/CMakeLists.txt
 rename examples/{helmholtzFlow/helmholtzDuctHeteroX => helmholtz2d/waveguide/longitudinal}/main.cpp (99%)
 rename examples/{helmholtzFlow/helmholtzDuctHeteroX/main_loop_length.cpp => helmholtz2d/waveguide/longitudinal/main_length.cpp} (99%)
 rename examples/{helmholtzFlow/helmholtzDuctHeteroX/main_loop_omega.cpp => helmholtz2d/waveguide/longitudinal/main_omega.cpp} (99%)
 create mode 100644 examples/helmholtz2d/waveguide/transverse/CMakeLists.txt
 rename examples/{helmholtzFlow/helmholtzDuctHeteroY => helmholtz2d/waveguide/transverse}/main.cpp (100%)
 rename examples/{helmholtzFlow/helmholtzDuctHeteroY/main_loop_length.cpp => helmholtz2d/waveguide/transverse/main_length.cpp} (100%)
 rename examples/{helmholtzFlow/helmholtzDuctHeteroY/main_loop_omega.cpp => helmholtz2d/waveguide/transverse/main_omega.cpp} (100%)
 create mode 100644 examples/helmholtzFlow/flowPastCylinder/CMakeLists.txt
 rename examples/helmholtzFlow/{helmholtzflow/FlowPastCylinder => flowPastCylinder}/main.cpp (100%)
 create mode 100644 examples/helmholtzFlow/freefield/CMakeLists.txt
 create mode 100644 examples/helmholtzFlow/freefield/main.cpp
 rename examples/helmholtzFlow/{helmholtzflow/FreeField => freefield}/runTests_circle.sh (72%)
 rename examples/helmholtzFlow/{helmholtzflow/FreeField => freefield}/runTests_square.sh (63%)
 delete mode 100644 examples/helmholtzFlow/helmholtzDuctHeteroX/CMakeLists.txt
 delete mode 100644 examples/helmholtzFlow/helmholtzDuctHeteroY/CMakeLists.txt
 delete mode 100644 examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/CMakeLists.txt
 delete mode 100644 examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/CMakeLists.txt
 create mode 100644 examples/helmholtzFlow/waveguide/Readme.md
 create mode 100644 examples/helmholtzFlow/waveguide/nonUniform/CMakeLists.txt
 rename examples/helmholtzFlow/{helmholtzflow/Duct_HeteroX => waveguide/nonUniform}/main.cpp (99%)
 rename examples/helmholtzFlow/{helmholtzflow/Duct_HeteroX/main_loop_omega.cpp => waveguide/nonUniform/main_omega.cpp} (100%)
 create mode 100644 examples/helmholtzFlow/waveguide/uniform/CMakeLists.txt
 create mode 100644 examples/helmholtzFlow/waveguide/uniform/main.cpp
 rename examples/helmholtzFlow/{helmholtzflow/Duct_ABCs => waveguide/uniform}/main_omega.cpp (100%)

diff --git a/contrib/bessel/CMakeLists.txt b/contrib/bessel/CMakeLists.txt
index 357acc31..898de89d 100644
--- a/contrib/bessel/CMakeLists.txt
+++ b/contrib/bessel/CMakeLists.txt
@@ -12,4 +12,4 @@ if(ENABLE_BESSEL)
 endif(ENABLE_BESSEL)
 
 file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
-append_gmshfem_contrib(Numeric "${SRC};${HDR}")
+append_gmshfem_contrib(bessel "${SRC};${HDR}")
diff --git a/examples/helmholtz2d/CMakeLists.txt b/examples/helmholtz2d/diskScattering/CMakeLists.txt
similarity index 100%
rename from examples/helmholtz2d/CMakeLists.txt
rename to examples/helmholtz2d/diskScattering/CMakeLists.txt
diff --git a/examples/helmholtz2d/README.md b/examples/helmholtz2d/diskScattering/README.md
similarity index 100%
rename from examples/helmholtz2d/README.md
rename to examples/helmholtz2d/diskScattering/README.md
diff --git a/examples/helmholtz2d/main.cpp b/examples/helmholtz2d/diskScattering/main.cpp
old mode 100755
new mode 100644
similarity index 100%
rename from examples/helmholtz2d/main.cpp
rename to examples/helmholtz2d/diskScattering/main.cpp
diff --git a/examples/helmholtz2d/waveguide/Readme.md b/examples/helmholtz2d/waveguide/Readme.md
new file mode 100644
index 00000000..67cc7e98
--- /dev/null
+++ b/examples/helmholtz2d/waveguide/Readme.md
@@ -0,0 +1,99 @@
+# Assessment of Local Absorbing Boundary Conditions for Heterogeneous Helmholtz Problems
+> Marchner P.
+
+Dependency(ies):
+* GmshFEM : v 1.0.0
+* Gmsh : v 4.8.2
+
+## Problem description
+
+This code solves a 2d straight waveguide Helmholtz problem 
+$$\rho_0^{-1}\partial_x(\rho_0\partial_x) u + \rho_0^{-1}\partial_y(\rho_0\partial_y) u + \omega^2 c_0^{-2} u = 0, \quad \text{in} \; \Omega=[0, L]\times[0, H]$$
+with a either 
+* a longitudinal, linear speed of sound profile $$c_0(x)=ax+b, \, a>0$$ (`longitudinal/`), 
+* a transverse, Gaussian $$c_0(y) ,\, \rho_0(y)$$ variation of the speed of sound and density (`transverse/`).
+
+High-order finite elements are used to discretize the weak formulations.
+A given mode of the form $$g=\cos(n \pi y / H), n\in\mathbb{N}$$ is imposed as input boundary condition at $$x=0$$.  
+Homogeneous Dirichlet or Neumann boundary conditions are available for the upper and lower walls.
+
+Different ABCs are tested on the output boundary.
+* In the folder `longitudinal/`, the relative $$L^2$$-error is computed thanks to an analytic solution.
+* In the folder `transverse/`, the relative $$L^2$$-error is computed thanks to a PML reference solution. 
+For the PML reference, Bermudez' function is used with the parameter $$\sigma_0=40$$ and 40 layers in order to ensure the attenuation of all modes including grazing modes.
+
+## Installation
+
+For the examples in the folder `longitudinal`, GmshFEM must be compiled with the Complex Bessel library
+
+```
+cmake -DENABLE_BESSEL=1 ..
+```
+
+Once GmshFEM is installed, go to the folder and run
+
+```
+  mkdir build && cd build
+  cmake ..
+  make
+```
+
+## Usage
+
+Simply run:
+```
+  ./example [PARAM]
+```
+with `[PARAM]`:
+* `-h [x]` where `[x]` is the mesh size.
+* `-L [x]` is the waveguide length and ABC position.
+* `-H [x]` is the waveguide height.
+* `-FEMorder [x]` is the order of the finite element hierarchical basis.
+* `-problem [x]` could be `hard` or `soft` and corresponds respectively to homogeneous Neumann and Dirichlet boundary conditions for the upper and lower duct walls.
+* `-omega [x]` is the angular frequency.
+* `-mode [x]` is the mode number for the input boundary condition.
+* `-abcName [x]` is the ABC imposed on the right boundary.
+* `-padeOrder [x]` is the Pade order if `abcName = Pade`. 
+* `-angle [x]` is the rotation angle of the branch-cut if `abcName = Pade`.
+* `-runcase [x]` selects a specific configuration for the results reproductibility.
+
+Each subfolder contains
+* the main script (`main.cpp`)
+* a loop for testing the ABCs over a given frequency range (`main_omega.cpp`)
+* a loop for varying the position of the ABC (`main_length.cpp`)
+To run one of the loops just change the name in the `CMakeLists.txt` and recompile
+
+## Results reproductibility
+##### Folder `longitudinal/` reproduce the results from Section 3
+Pressure maps from Figure 2 can be obtained by compiling `main.cpp` and running for example
+```
+  ./example -omega 20
+```
+The solution `v_exact.msh` can be open in Gmsh.
+
+For Figure 3 use `main_omega.cpp` and run
+```
+  ./example -runcase [i]
+```
+with i=1 for Taylor-based ABCs or i=2 for Padé-based ABCs.
+
+For Figure 4 use `main_length.cpp`
+```
+  ./example -runcase [i]
+```
+with i={1,2,3} to obtain the data from Figures 4a, 4b and 4c respectively.
+
+##### Folder `transverse/` reproduce the results from Section 5
+For Figures 13a, 13b and 13c, use respectively `main.cpp`, `main_omega.cpp` and `main_length.cpp` with
+```
+./example -runcase 1
+```
+
+Figure 15a and 15b, use respectively `main.cpp` and `main_omega.cpp` with
+```
+./example -runcase 2
+```
+
+## Reference
+
+// Preprint
diff --git a/examples/helmholtz2d/waveguide/longitudinal/CMakeLists.txt b/examples/helmholtz2d/waveguide/longitudinal/CMakeLists.txt
new file mode 100644
index 00000000..6f8caaab
--- /dev/null
+++ b/examples/helmholtz2d/waveguide/longitudinal/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(example CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake)
+
+add_executable(example main.cpp ${EXTRA_INCS})
+target_link_libraries(example ${EXTRA_LIBS})
diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroX/main.cpp b/examples/helmholtz2d/waveguide/longitudinal/main.cpp
similarity index 99%
rename from examples/helmholtzFlow/helmholtzDuctHeteroX/main.cpp
rename to examples/helmholtz2d/waveguide/longitudinal/main.cpp
index fda2673c..0661c47d 100644
--- a/examples/helmholtzFlow/helmholtzDuctHeteroX/main.cpp
+++ b/examples/helmholtz2d/waveguide/longitudinal/main.cpp
@@ -76,7 +76,6 @@ int main(int argc, char **argv)
   int FEMorder = 6;
   gmshFem.userDefinedParameter(FEMorder, "FEMorder");
   std::string gauss = "Gauss"+std::to_string(2*FEMorder+1);
-  gmshFem.userDefinedParameter(gauss, "gauss");
 
   // meshing
   double L = 1; // duct length
diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_length.cpp b/examples/helmholtz2d/waveguide/longitudinal/main_length.cpp
similarity index 99%
rename from examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_length.cpp
rename to examples/helmholtz2d/waveguide/longitudinal/main_length.cpp
index c953106e..2bce8fb8 100644
--- a/examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_length.cpp
+++ b/examples/helmholtz2d/waveguide/longitudinal/main_length.cpp
@@ -80,7 +80,6 @@ int main(int argc, char **argv)
   int FEMorder = 6;
   gmshFem.userDefinedParameter(FEMorder, "FEMorder");
   std::string gauss = "Gauss"+std::to_string(2*FEMorder+1);
-  gmshFem.userDefinedParameter(gauss, "gauss");
 
   // meshing
   double H = 0.5; // duct height
diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_omega.cpp b/examples/helmholtz2d/waveguide/longitudinal/main_omega.cpp
similarity index 99%
rename from examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_omega.cpp
rename to examples/helmholtz2d/waveguide/longitudinal/main_omega.cpp
index a2f302ac..73421a43 100644
--- a/examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_omega.cpp
+++ b/examples/helmholtz2d/waveguide/longitudinal/main_omega.cpp
@@ -79,7 +79,6 @@ int main(int argc, char **argv)
   int FEMorder = 6;
   gmshFem.userDefinedParameter(FEMorder, "FEMorder");
   std::string gauss = "Gauss"+std::to_string(2*FEMorder+1);
-  gmshFem.userDefinedParameter(gauss, "gauss");
 
   // meshing
   double L = 1; // duct length
diff --git a/examples/helmholtz2d/waveguide/transverse/CMakeLists.txt b/examples/helmholtz2d/waveguide/transverse/CMakeLists.txt
new file mode 100644
index 00000000..6f8caaab
--- /dev/null
+++ b/examples/helmholtz2d/waveguide/transverse/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(example CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake)
+
+add_executable(example main.cpp ${EXTRA_INCS})
+target_link_libraries(example ${EXTRA_LIBS})
diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroY/main.cpp b/examples/helmholtz2d/waveguide/transverse/main.cpp
similarity index 100%
rename from examples/helmholtzFlow/helmholtzDuctHeteroY/main.cpp
rename to examples/helmholtz2d/waveguide/transverse/main.cpp
diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroY/main_loop_length.cpp b/examples/helmholtz2d/waveguide/transverse/main_length.cpp
similarity index 100%
rename from examples/helmholtzFlow/helmholtzDuctHeteroY/main_loop_length.cpp
rename to examples/helmholtz2d/waveguide/transverse/main_length.cpp
diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroY/main_loop_omega.cpp b/examples/helmholtz2d/waveguide/transverse/main_omega.cpp
similarity index 100%
rename from examples/helmholtzFlow/helmholtzDuctHeteroY/main_loop_omega.cpp
rename to examples/helmholtz2d/waveguide/transverse/main_omega.cpp
diff --git a/examples/helmholtzFlow/flowPastCylinder/CMakeLists.txt b/examples/helmholtzFlow/flowPastCylinder/CMakeLists.txt
new file mode 100644
index 00000000..6f8caaab
--- /dev/null
+++ b/examples/helmholtzFlow/flowPastCylinder/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(example CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake)
+
+add_executable(example main.cpp ${EXTRA_INCS})
+target_link_libraries(example ${EXTRA_LIBS})
diff --git a/examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/main.cpp b/examples/helmholtzFlow/flowPastCylinder/main.cpp
similarity index 100%
rename from examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/main.cpp
rename to examples/helmholtzFlow/flowPastCylinder/main.cpp
diff --git a/examples/helmholtzFlow/freefield/CMakeLists.txt b/examples/helmholtzFlow/freefield/CMakeLists.txt
new file mode 100644
index 00000000..93e63bab
--- /dev/null
+++ b/examples/helmholtzFlow/freefield/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(example CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../example.cmake)
+
+add_executable(example main.cpp ${EXTRA_INCS})
+target_link_libraries(example ${EXTRA_LIBS})
diff --git a/examples/helmholtzFlow/freefield/main.cpp b/examples/helmholtzFlow/freefield/main.cpp
new file mode 100644
index 00000000..a63b33ac
--- /dev/null
+++ b/examples/helmholtzFlow/freefield/main.cpp
@@ -0,0 +1,564 @@
+#include "GmshFem.h"
+#include "Formulation.h"
+#include "gmsh.h"
+#include "AnalyticalFunction.h"
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::problem;
+using namespace gmshfem::domain;
+using namespace gmshfem::field;
+using namespace gmshfem::function;
+using namespace gmshfem::analytics;
+using namespace gmshfem::post;
+using namespace gmshfem::equation;
+
+struct Edge
+{
+  Domain gamma;
+  Domain corner[2];
+};
+
+struct HelmholtzDomain
+{
+  Domain omega;
+  Domain gammaPhy;
+  Domain gammaPml;
+  Domain omegaPml;
+  Domain source;
+  Domain corners;
+  std::vector< Domain > corner;
+  std::vector< Edge > edge;
+};
+
+//***************************************
+// GEOMETRY
+//***************************************
+
+void meshCircle(GmshFem &gmshFem, HelmholtzDomain &domains, const double R, const double Rpml, const double xs, const double ys, const double lc, bool withPml)
+{
+  gmsh::model::add("geometry");
+
+  gmsh::model::geo::addPoint(0., 0., 0., lc, 1); // Center
+  gmsh::model::geo::addPoint(xs, ys, 0., lc, 100); // Source point
+  // physical circle
+  gmsh::model::geo::addPoint(R, 0., 0., lc, 2);
+  gmsh::model::geo::addPoint(0, R, 0., lc, 3);
+  gmsh::model::geo::addPoint(-R, 0., 0., lc, 4);
+  gmsh::model::geo::addPoint(0, -R, 0., lc, 5);
+
+  gmsh::model::geo::addCircleArc(2, 1, 3, 1);
+  gmsh::model::geo::addCircleArc(3, 1, 4, 2);
+  gmsh::model::geo::addCircleArc(4, 1, 5, 3);
+  gmsh::model::geo::addCircleArc(5, 1, 2, 4);
+  gmsh::model::geo::addCurveLoop({1,2,3,4}, 1);
+  // surface
+  gmsh::model::geo::addPlaneSurface({1}, 1);
+
+  if(withPml) {
+    gmsh::model::geo::addPoint(Rpml, 0., 0., lc, 6);
+    gmsh::model::geo::addPoint(0, Rpml, 0., lc, 7);
+    gmsh::model::geo::addPoint(-Rpml, 0., 0., lc, 8);
+    gmsh::model::geo::addPoint(0, -Rpml, 0., lc, 9);
+
+    gmsh::model::geo::addCircleArc(6, 1, 7, 5);
+    gmsh::model::geo::addCircleArc(7, 1, 8, 6);
+    gmsh::model::geo::addCircleArc(8, 1, 9, 7);
+    gmsh::model::geo::addCircleArc(9, 1, 6, 8);
+    gmsh::model::geo::addCurveLoop({5,6,7,8}, 2);
+    gmsh::model::geo::addPlaneSurface({1,2}, 2);
+  }
+  gmsh::model::geo::synchronize();
+  gmsh::model::mesh::embed(0,{100},2,1);
+
+  // physicals
+  gmsh::model::addPhysicalGroup(2, {1}, 1);
+  gmsh::model::setPhysicalName(2, 1, "omega");
+  gmsh::model::addPhysicalGroup(1, {1,2,3,4}, 3);
+  gmsh::model::setPhysicalName(1, 3, "gammaPhy");
+  gmsh::model::addPhysicalGroup(0, {100}, 100);
+  gmsh::model::setPhysicalName(0, 100, "source");
+  if(withPml) {
+    gmsh::model::addPhysicalGroup(2, {2}, 2);
+    gmsh::model::setPhysicalName(2, 2, "omegaPml");
+    gmsh::model::addPhysicalGroup(1, {5,6,7,8}, 4);
+    gmsh::model::setPhysicalName(1, 4, "gammaPml");
+  }
+
+  // generate mesh
+  gmsh::model::geo::synchronize();
+  gmsh::model::mesh::generate(2);
+  // gmsh::model::mesh::setOrder(2); // T6 elem
+  // gmsh::write("m.msh");
+  
+  // Domains
+  domains.omega = Domain(2, 1);
+  domains.gammaPhy = Domain(1, 3);
+  domains.source = Domain(0, 100);
+
+  if(withPml) {
+    domains.omegaPml = Domain(2, 2);
+    domains.gammaPml = Domain(1, 4);
+  }
+}
+
+void meshSquare(GmshFem &gmshFem, HelmholtzDomain &domains, const double L, const double Lpml, const double xs, const double ys, const double lc, bool withPml)
+{
+  gmsh::model::add("geometry");
+
+  gmsh::model::geo::addPoint(0., 0., 0., lc, 1); // Center
+  gmsh::model::geo::addPoint(xs, ys, 0., lc, 100); // Source point
+
+  // polygone
+  gmsh::model::geo::addPoint(L, L, 0., lc, 2);
+  gmsh::model::geo::addPoint(-L, L, 0., lc, 3);
+  gmsh::model::geo::addPoint(-L, -L, 0., lc, 4);
+  gmsh::model::geo::addPoint(L, -L, 0., lc, 5);
+
+  gmsh::model::geo::addLine(2, 3, 1);
+  gmsh::model::geo::addLine(3, 4, 2);
+  gmsh::model::geo::addLine(4, 5, 3);
+  gmsh::model::geo::addLine(5, 2, 4);
+
+  gmsh::model::geo::addCurveLoop({1,2,3,4}, 1);
+  gmsh::model::geo::addPlaneSurface({1}, 1);
+
+  if(withPml) {
+    gmsh::model::geo::addPoint(Lpml, Lpml, 0., lc, 6);
+    gmsh::model::geo::addPoint(-Lpml, Lpml, 0., lc, 7);
+    gmsh::model::geo::addPoint(-Lpml, -Lpml, 0., lc, 8);
+    gmsh::model::geo::addPoint(Lpml, -Lpml, 0., lc, 9);
+
+    gmsh::model::geo::addLine(6, 7, 5);
+    gmsh::model::geo::addLine(7, 8, 6);
+    gmsh::model::geo::addLine(8, 9, 7);
+    gmsh::model::geo::addLine(9, 6, 8);
+
+    gmsh::model::geo::addCurveLoop({5,6,7,8}, 2);
+    gmsh::model::geo::addPlaneSurface({1,2}, 2);
+  }
+  gmsh::model::geo::synchronize();
+  gmsh::model::mesh::embed(0,{100},2,1);
+  // physicals
+  gmsh::model::addPhysicalGroup(2, {1}, 1);
+  gmsh::model::setPhysicalName(2, 1, "omega");
+  gmsh::model::addPhysicalGroup(1, {1,2,3,4}, 3);
+  gmsh::model::setPhysicalName(1, 3, "gammaPhy");
+  gmsh::model::addPhysicalGroup(0, {100}, 100);
+  gmsh::model::setPhysicalName(0, 100, "source");
+  // Physical edges
+  gmsh::model::addPhysicalGroup(1, {1}, 20);
+  gmsh::model::setPhysicalName(1, 20, "gammaPhy_top");
+  gmsh::model::addPhysicalGroup(1, {2}, 30);
+  gmsh::model::setPhysicalName(1, 30, "gammaPhy_left");
+  gmsh::model::addPhysicalGroup(1, {3}, 40);
+  gmsh::model::setPhysicalName(1, 40, "gammaPhy_bottom");
+  gmsh::model::addPhysicalGroup(1, {4}, 50);
+  gmsh::model::setPhysicalName(1, 50, "gammaPhy_right");
+  // Physical domain corners
+  gmsh::model::addPhysicalGroup(0, {2}, 2);
+  gmsh::model::setPhysicalName(0, 2, "UpRightCorner");
+  gmsh::model::addPhysicalGroup(0, {3}, 3);
+  gmsh::model::setPhysicalName(0, 3, "UpLeftCorner");
+  gmsh::model::addPhysicalGroup(0, {4}, 4);
+  gmsh::model::setPhysicalName(0, 4, "BottomLeftCorner");
+  gmsh::model::addPhysicalGroup(0, {5}, 5);
+  gmsh::model::setPhysicalName(0, 5, "BottomRightCorner");
+  if(withPml) {
+    gmsh::model::addPhysicalGroup(2, {2}, 2);
+    gmsh::model::setPhysicalName(2, 2, "omegaPml");
+    gmsh::model::addPhysicalGroup(1, {5,6,7,8}, 4);
+    gmsh::model::setPhysicalName(1, 4, "gammaPml");
+  }
+
+  // generate mesh
+  gmsh::model::geo::synchronize();
+  gmsh::model::mesh::generate(2);
+  //gmsh::model::mesh::setOrder(2); // T6 elem
+  // gmsh::write("m.msh");
+  
+  // Domains
+  domains.omega = Domain(2, 1);
+  domains.gammaPhy = Domain(1, 3);
+  domains.source = Domain(0, 100);
+
+  if(withPml) {
+    domains.omegaPml = Domain(2, 2);
+    domains.gammaPml = Domain(1, 4);
+  }
+  domains.edge.resize(4);
+  domains.corner.resize(4);
+  for(int i = 0; i < 4; ++i) {
+    domains.edge[i].gamma = Domain(1, 10*(i + 2) );
+    domains.edge[i].corner[0] = Domain(0, i + 2 );
+    domains.edge[i].corner[1] = Domain(0, (i + 1) % 4 + 2 );
+    
+    domains.corner[i] = Domain(0, i + 2);
+    
+    domains.corners |= domains.edge[i].corner[0];
+    domains.gammaPhy |= domains.edge[i].gamma;
+  }
+}
+
+//***************************************
+// FORMULATION
+//***************************************
+
+int main(int argc, char **argv)
+{
+  GmshFem gmshFem(argc, argv);
+  
+  // geometry parameters
+  int geometry = 0; // 0-circle, 1-square
+  gmshFem.userDefinedParameter(geometry, "geometry");
+  double R = 1.0; // Circle radius or square of size [-R,R] x [-R,R]
+  gmshFem.userDefinedParameter(R, "R");
+  // Point source location
+  double xs = 0.3;
+  double ys = 0.2;
+  gmshFem.userDefinedParameter(xs, "xs");
+  gmshFem.userDefinedParameter(ys, "ys");
+  msg::info << "Point source located at xs=" << xs << ", ys=" << ys << msg::endl;
+  
+  // physical parameters
+  const double pi = 3.14159265358979323846264338327950288;
+  const double k = 3 * pi; // free field wavenumber
+  const std::complex< double > im(0., 1.);
+  double M = 0.; // Mach number
+  gmshFem.userDefinedParameter(M, "M");
+  double theta = pi/4; // mean flow orientation
+  gmshFem.userDefinedParameter(theta, "theta");
+  msg::info << "Mach number " << M << " oriented at theta=" << theta << " rad" << msg::endl;
+  
+  // numerical parameters
+  int order = 4; // FEM shape function order
+  gmshFem.userDefinedParameter(order, "order");
+  std::string gauss = "Gauss10"; // increase Nr of Gauss points when order > 6
+  double lc = 0.03; // mesh size
+  gmshFem.userDefinedParameter(lc, "lc");
+  
+  // Choose ABC - exterior boundary condition
+  std::string abcName = "Pade"; // available choices: ABC-0, ABC-2, Pade, pml
+  gmshFem.userDefinedParameter(abcName, "abcName");
+  int Npml;
+  double Rpml;
+  if (abcName == "pml") { // PML parameters
+    Npml = 4; // number of layers
+    gmshFem.userDefinedParameter(Npml, "Npml");
+    Rpml = R+Npml*lc; // extended domain
+  }
+
+  HelmholtzDomain domains;
+
+  switch (geometry) {
+  case 0:
+    meshCircle(gmshFem, domains, R, Rpml, xs, ys, lc, abcName == "pml");
+    break;
+  case 1:
+    meshSquare(gmshFem, domains, R, Rpml, xs, ys, lc, abcName == "pml");
+    break;
+  default:
+    msg::error << "geometry not available ! " << msg::endl;
+    exit(0);
+    break;
+  }
+  
+  float pointsByWl = (2*pi*order*(1-abs(M)))/(lc*k);
+  msg::info << " - FEM basis order = " << order << "" << msg::endl;
+  msg::info << " - Smallest number of dofs by wavelength = " << pointsByWl << "" << msg::endl;
+  if (pointsByWl < 6) {
+    msg::warning << " - Less than 6 points per wavelength ! " << msg::endl;
+  }
+
+  // Mach number projections 
+  double Mx = M*std::cos(theta);
+  double My = M*std::sin(theta);
+  // Mach-velocity vector
+  VectorFunction< std::complex< double > > MM = vector< std::complex< double > > (Mx,My,0.);
+  // Normal and tangential projections 
+  ScalarFunction< std::complex< double > > Mn = vector< std::complex< double > > (Mx,My,0.)*normal< std::complex< double > >();
+  ScalarFunction< std::complex< double > > Mt = vector< std::complex< double > > (Mx,My,0.)*tangent< std::complex< double > >();
+  // Parameters for Inverse Lorentz transformation
+  double beta = sqrt(1-M*M); // Jacobian
+  double Alphax = 1 + Mx*Mx/(beta*(1+beta));
+  double Alphay = 1 + My*My/(beta*(1+beta));
+  double K = Mx*My/(beta*(1+beta));
+  // Tensor of the Inverse Lorentz transformation
+  TensorFunction< std::complex< double > > Linv = tensor< std::complex <double > > (beta*Alphay, -beta*K, 0., -beta*K, beta*Alphax,0.,0.,0.,0.);
+
+  std::vector< FieldInterface< std::complex< double > > * > fieldBucket;
+  Formulation< std::complex< double > > formulation("helmholtzflow");
+  Field< std::complex< double >, Form::Form0 > u("u", domains.omega | domains.omegaPml | domains.gammaPhy | domains.gammaPml | domains.corners | domains.source, FunctionSpaceTypeForm0::HierarchicalH1, order);
+
+  // Define analytic solution
+  Function< std::complex< double >, Degree::Degree0 > *solution = nullptr;
+  solution = new AnalyticalFunction< helmholtz2D::MonopoleFreeField< std::complex< double > > >(k, M, theta, 1., xs, ys, 0., 0.);
+  solution->activateMemorization();
+
+  // convected Helmholz weak form
+  formulation.galerkin(vector< std::complex< double > >(1-Mx*Mx, -Mx*My, 0.) * grad(dof(u)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(u)), domains.omega, gauss);
+  formulation.galerkin(vector< std::complex< double > >(-Mx*My, 1-My*My, 0.) * grad(dof(u)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omega, gauss);
+  formulation.galerkin(- k * k * dof(u), tf(u), domains.omega, gauss);
+  formulation.galerkin(dof(u), vector< std::complex< double > >(-im*k*Mx,-im*k*My, 0.) * grad(tf(u)), domains.omega, gauss);
+  formulation.galerkin(vector< std::complex< double > >(im*k*Mx, im*k*My, 0.) * grad(dof(u)), tf(u), domains.omega, gauss);
+
+  // Source point
+  formulation.galerkin(-1.,tf(u),domains.source,gauss);
+
+  // Exact DtN: beta_n**2 * kx = - Mn * (k0 - Mt*ky) + sqrt( k0**2 -2*k0*Mt*ky - (1-abs(M))^2 ky**2 ) 
+  if(abcName == "ABC-0") { // zeroth order Taylor approximation
+    // boundary contributions
+    formulation.galerkin(im*k*Mn * dof(u), tf(u), domains.gammaPhy, gauss);
+    formulation.galerkin(Mn*Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss);
+    // ABC
+    msg::info << "Use zeroth order ABC" << msg::endl;
+    formulation.galerkin(im*k*(1-Mn) * dof(u), tf(u), domains.gammaPhy, gauss);
+    formulation.galerkin(-Mn*Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss);
+  }
+  else if (abcName == "ABC-2") {
+    // boundary contributions
+    formulation.galerkin(im*k*Mn * dof(u), tf(u), domains.gammaPhy, gauss);
+    formulation.galerkin(Mn*Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss);
+    // ABC
+    msg::info << "Use second order ABC" << msg::endl;
+    // zeroth order contribution
+    formulation.galerkin(im*k*(1-Mn) * dof(u), tf(u), domains.gammaPhy, gauss);
+    // first order contribution
+    formulation.galerkin( (1-Mn) * Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss);
+    // second order contribution
+    formulation.galerkin(- im * (beta*beta) /(2*k) * grad(dof(u)), grad(tf(u)), domains.gammaPhy, gauss);
+  }
+  else if (abcName == "Pade") {
+    int padeOrder = 2;
+    gmshFem.userDefinedParameter(padeOrder, "padeOrder");
+    double angle = -pi/4.;
+    gmshFem.userDefinedParameter(angle, "angle");
+    msg::info << "Use Pade ABC of order " << padeOrder << " with angle " << angle << " rad"<< msg::endl;
+    
+    const double Np = 2. * padeOrder + 1.;
+    const std::complex< double > exp1 = std::complex<double>(std::cos(angle),std::sin(angle));
+    const std::complex< double > exp2 = std::complex<double>(std::cos(angle/2.),std::sin(angle/2.));
+    const std::complex< double > coef = 2./Np;
+    std::vector< std::complex< double > > c(padeOrder, 0.);
+    for(int i = 0; i < padeOrder; ++i) {
+      c[i] = std::tan((i + 1) * pi / Np);
+      c[i] *= c[i];
+    }
+    
+    if(geometry == 0) {
+      // define the auxiliary fields
+      std::vector< Field< std::complex< double >, Form::Form0 >* > phi;
+      for(int i = 0; i < padeOrder; ++i) {
+        phi.push_back(new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(i), domains.gammaPhy, FunctionSpaceTypeForm0::HierarchicalH1, order));
+        fieldBucket.push_back(phi.back());
+      }
+        
+      // write the augmented weak form - approximation of the square-root
+      formulation.galerkin( im * k * exp2 * dof(u), tf(u), domains.gammaPhy, gauss);
+      for(int i = 0; i < padeOrder; ++i) {
+        // boundary integral terms relating the auxiliary fields
+        formulation.galerkin( im * k * exp2 * coef * c[i] * dof(*phi[i]), tf(u), domains.gammaPhy, gauss);
+        formulation.galerkin( im * k * exp2 * coef * c[i] * dof(u), tf(u), domains.gammaPhy, gauss);
+
+        // coupling of the auxiliary equations
+        formulation.galerkin(- (beta*beta) * grad(dof(*phi[i])), grad(tf(*phi[i])), domains.gammaPhy, gauss);
+        formulation.galerkin(- 2. * im * k * Mt * tangent< std::complex< double > >() * grad(dof(*phi[i])), tf(*phi[i]), domains.gammaPhy, gauss);
+        formulation.galerkin( (k*k) * (exp1 * c[i] + 1.) * dof(*phi[i]), tf(*phi[i]), domains.gammaPhy, gauss);
+        formulation.galerkin( (k*k) * exp1 * (c[i] + 1.) * dof(u), tf(*phi[i]), domains.gammaPhy, gauss);
+      }
+    }
+    else if(geometry == 1) {
+      // define the auxilary fields
+      std::vector< Field< std::complex< double >, Form::Form0 > * > phi[4];
+      for(int x = 0; x < 4; ++x) {
+        for(int i = 0; i < padeOrder; ++i) {
+          phi[x].push_back(new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(i) + "^" + std::to_string(x), domains.edge[x].gamma | domains.edge[x].corner[0] | domains.edge[x].corner[1], FunctionSpaceTypeForm0::HierarchicalH1,order));
+          fieldBucket.push_back(phi[x].back());
+        }
+      }
+      
+      formulation.galerkin( im * k * exp2 * dof(u), tf(u), domains.gammaPhy, gauss);
+      for(int x = 0; x < 4; ++x) {
+        for(int i = 0; i < padeOrder; ++i) {
+          // Boundary conditions
+          formulation.galerkin( im * k * exp2 * coef * c[i] * dof(*phi[x][i]), tf(u), domains.edge[x].gamma, gauss);
+          formulation.galerkin( im * k * exp2 * coef * c[i] * dof(u), tf(u), domains.edge[x].gamma, gauss);
+          
+          // Auxiliary equations
+          formulation.galerkin(- (beta*beta) * grad(dof(*phi[x][i])), grad(tf(*phi[x][i])), domains.edge[x].gamma, gauss);
+          formulation.galerkin(- 2. * im * k * Mt * tangent< std::complex< double > >() * grad(dof(*phi[x][i])), tf(*phi[x][i]), domains.edge[x].gamma, gauss);
+          formulation.galerkin( (k*k) * (exp1 * c[i] + 1.) * dof(*phi[x][i]), tf(*phi[x][i]), domains.edge[x].gamma, gauss);
+          formulation.galerkin( (k*k) * exp1 * (c[i] + 1.) * dof(u), tf(*phi[x][i]), domains.edge[x].gamma, gauss);
+        }
+      }
+      
+      bool withCorner = true;
+      gmshFem.userDefinedParameter(withCorner, "withCorner");
+      std::string CornerStrategy = "Sommerfeld"; // Sommerfeld or HABC 
+      gmshFem.userDefinedParameter(CornerStrategy, "CornerStrategy");
+      
+      if(withCorner) {
+        msg::info << "with corner treatment" << msg::endl;
+        msg::info << "use " << CornerStrategy << " condition at corners" << msg::endl;
+        // define the corner fields
+        std::vector< std::vector< std::vector< Field< std::complex< double >, Form::Form0 > * > > > corner_phi(4);
+        for(int k = 0; k < 4; ++k) {
+          corner_phi[k].resize(padeOrder);
+          for(int i = 0; i < padeOrder; ++i) {
+            for(int j = 0; j < padeOrder; ++j) {
+              corner_phi[k][i].push_back(new Field< std::complex< double >, Form::Form0 >("corner_phi_(" + std::to_string(i) + ", " + std::to_string(j) + ")^(" + std::to_string(k) + ")", domains.corner[k], FunctionSpaceTypeForm0::HierarchicalH1,order));
+              fieldBucket.push_back(corner_phi[k][i].back());
+            }
+          }
+        }
+        
+        for(int ci = 0; ci < 4; ++ci) {
+          int x = ci;
+          int y = (ci - 1 < 0 ? 4 - 1 : ci - 1);
+          for(int i = 0; i < padeOrder; ++i) {
+            if (CornerStrategy == "Sommerfeld"){
+              // strategy 1 - Sommerfeld condition at corners
+              formulation.galerkin(- im * k * dof(*phi[x][i]), tf(*phi[x][i]), domains.corner[ci], gauss);
+              formulation.galerkin(- im * k * dof(*phi[y][i]), tf(*phi[y][i]), domains.corner[ci], gauss);
+            }
+
+            else if (CornerStrategy == "HABC") {
+              // strategy 2 - HABC at corners
+              formulation.galerkin(- im * k * exp2 * dof(*phi[x][i]), tf(*phi[x][i]), domains.corner[ci], gauss);
+              formulation.galerkin(- im * k * exp2 * dof(*phi[y][i]), tf(*phi[y][i]), domains.corner[ci], gauss);
+              for(int j = 0; j < padeOrder; ++j) {
+                // Auxiliary equations for corner
+                formulation.galerkin( (exp1*c[i] + exp1*c[j] + 1.) * dof(*corner_phi[ci][i][j]), tf(*corner_phi[ci][i][j]), domains.corner[ci], gauss);
+                formulation.galerkin( exp1 * (c[j] + 1.) * dof(*phi[x][i]), tf(*corner_phi[ci][i][j]), domains.corner[ci], gauss);
+                formulation.galerkin( exp1 * (c[i] + 1.) * dof(*phi[y][j]), tf(*corner_phi[ci][i][j]), domains.corner[ci], gauss);
+                
+                // Corner conditions
+                formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*phi[x][i]), tf(*phi[x][i]), domains.corner[ci], gauss);
+                formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*phi[y][i]), tf(*phi[y][i]), domains.corner[ci], gauss);
+                formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*corner_phi[ci][i][j]), tf(*phi[x][i]), domains.corner[ci], gauss);
+                formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*corner_phi[ci][j][i]), tf(*phi[y][i]), domains.corner[ci], gauss);
+              }
+            }
+            else {
+              msg::error << "Corner strategy not available !" << msg::endl;
+              exit(0);
+            }
+          }
+        }
+      }
+    }
+  }
+  else if (abcName == "pml") {
+    msg::info << "Use a PML with " << Npml << " layers" << msg::endl;
+    double Sigma0 = beta;
+    ScalarFunction< std::complex< double > > det_J;
+    TensorFunction< std::complex< double > > J_PML_inv_T;
+    const double Wpml = Rpml - R;
+    if (geometry == 0) {
+      msg::info << "stabilized PML - circular domain" << msg::endl;
+      msg::info << "use unbounded absorbing function with sigma_0 = " << Sigma0 << msg::endl;
+      ScalarFunction< std::complex< double > > cosT = x< std::complex< double > >() / r2d< std::complex< double > >();
+      ScalarFunction< std::complex< double > > sinT = y< std::complex< double > >() / r2d< std::complex< double > >();
+      ScalarFunction< std::complex< double > > dampingProfileR = Sigma0 / (Wpml - (r2d< std::complex< double > >() - R));
+      ScalarFunction< std::complex< double > > dampingProfileInt = -Sigma0*ln((Wpml - (r2d< std::complex< double > >() - R)) / Wpml);
+      ScalarFunction< std::complex< double > > gamma = 1. - im * dampingProfileR / k;
+      ScalarFunction< std::complex< double > > gamma_hat = 1. - im * (1. / r2d< std::complex< double > >()) * dampingProfileInt / k;
+      det_J = gamma * gamma_hat;
+      J_PML_inv_T = tensor< std::complex <double > >(cosT/gamma,sinT/gamma,0.,-sinT/gamma_hat,cosT/gamma_hat,0.,0.,0.,0.);
+    }
+    else if (geometry == 1) {
+      msg::info << "stabilized PML - square domain" << msg::endl;
+      msg::info << "use unbounded absorbing function with sigma_0 = " << Sigma0 << msg::endl;
+      ScalarFunction< std::complex< double > > SigmaX = heaviside(abs(x< std::complex< double > >()) - R ) * Sigma0/(Rpml - abs(x< std::complex< double > >()) );
+      ScalarFunction< std::complex< double > > gammaX = 1-(im/k)*SigmaX;
+      ScalarFunction< std::complex< double > > SigmaY = heaviside(abs(y< std::complex< double > >()) - R ) * Sigma0/(Rpml - abs(y< std::complex< double > >()) );
+      ScalarFunction< std::complex< double > > gammaY = 1-(im/k)*SigmaY;
+      det_J = gammaX * gammaY;
+      J_PML_inv_T = tensor< std::complex <double > >(1./gammaX,0.,0.,0.,1./gammaY,0.,0.,0.,0.);
+
+      /*
+      msg::info << "classical PML - square domain - unstable ! " << msg::endl;
+      formulation.galerkin(vector< std::complex< double > >((gammaY/gammaX)*(1-Mx*Mx), -Mx*My, 0.) * grad(dof(u)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(u)), domains.omegaPml, gauss);
+      formulation.galerkin(vector< std::complex< double > >(-Mx*My, (gammaX/gammaY)*(1-My*My), 0.) * grad(dof(u)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omegaPml, gauss);
+      formulation.galerkin(dof(u), vector< std::complex< double > >(-gammaY*im*k*Mx,-gammaX*im*k*My, 0.) * grad(tf(u)), domains.omegaPml, gauss);
+      formulation.galerkin(vector< std::complex< double > >(gammaY*im*k*Mx, gammaX*im*k*My, 0.) * grad(dof(u)), tf(u), domains.omegaPml, gauss);
+      formulation.galerkin(- k*k* gammaX*gammaY * dof(u), tf(u), domains.omegaPml, gauss);
+      */
+
+      /*
+      msg::info << "Alternative stable PML - square domain only " << msg::endl;
+      ScalarFunction< std::complex< double > > Ax = (gammaY/gammaX)*(1-Mx*Mx);
+      std::complex< double > Cx =-im*k*Mx*( 1 - (My*My)/(1-Mx*Mx) ) / (beta*beta);
+      double Ka = -Mx*My/(1-Mx*Mx);
+      ScalarFunction< std::complex< double > > Ay = (gammaX/gammaY)*( 1 - (My*My)/(1-Mx*Mx) );
+      std::complex< double > Cy = -im*k*My/(beta*beta);
+
+      // modified bilinear term in dx
+      formulation.galerkin( Ax * vector< std::complex< double > >(1. , Ka, 0.) * grad(dof(u)), vector< std::complex< double > >(1., Ka, 0.) * grad(tf(u)), domains.omegaPml, gauss);
+      formulation.galerkin( Ax * Cx * dof(u), vector< std::complex< double > >(1., Ka, 0.) * grad(tf(u)), domains.omegaPml, gauss);
+      formulation.galerkin( Ax * (-Cx) * vector< std::complex< double > >(1. , Ka, 0.) * grad(dof(u)), tf(u), domains.omegaPml, gauss);
+      formulation.galerkin( Ax * Cx * (-Cx) *dof(u), tf(u), domains.omegaPml, gauss);
+      // modified bilinear term in dy
+      formulation.galerkin( Ay * vector< std::complex< double > >(0. , 1., 0.) * grad(dof(u)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omegaPml, gauss);
+      formulation.galerkin( Ay * Cy * dof(u), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omegaPml, gauss);
+      formulation.galerkin( Ay * vector< std::complex< double > >(0. , -Cy, 0.) * grad(dof(u)), tf(u), domains.omegaPml, gauss);
+      formulation.galerkin( Ay * Cy * (-Cy) *dof(u), tf(u), domains.omegaPml, gauss);
+      //
+      formulation.galerkin(- k*k/(beta*beta) * gammaX*gammaY * dof(u), tf(u), domains.omegaPml, gauss);
+      */
+    }
+    // build vector and matrix for the weak form
+    VectorFunction< std::complex< double > > J_PML_inv_T_M = J_PML_inv_T * MM;
+    TensorFunction< std::complex< double > > J_PML_Linv  = J_PML_inv_T * Linv;
+    // stabilized PML weak form - valid for a general domain
+    formulation.galerkin( det_J * J_PML_Linv*grad(dof(u)) , J_PML_Linv*grad(tf(u)) , domains.omegaPml, gauss);
+
+    formulation.galerkin(+ k*k/(beta*beta) * det_J * J_PML_inv_T_M * dof(u) , J_PML_inv_T_M * tf(u), domains.omegaPml, gauss);
+    formulation.galerkin(+ im*k/beta* det_J * J_PML_Linv * grad(dof(u)), J_PML_inv_T_M * tf(u), domains.omegaPml, gauss);
+
+    formulation.galerkin(- im*k/beta* det_J * J_PML_inv_T_M * dof(u), J_PML_Linv * grad(tf(u)), domains.omegaPml, gauss);
+    formulation.galerkin(- k*k/(beta*beta) * det_J * dof(u), tf(u), domains.omegaPml, gauss);
+  }
+  else {
+    msg::error << "ABC Type not available !" << msg::endl;
+    exit(0);
+  }
+
+  // solve
+  formulation.pre();
+  formulation.assemble();
+  formulation.solve();
+
+  // compute the projection of the analytic solution on the FEM basis - best FEM approximation
+  Formulation< std::complex< double > > projection("helmholtzflow");
+  Field< std::complex< double >, Form::Form0 > uP("uP", domains.omega, FunctionSpaceTypeForm0::HierarchicalH1, order);
+  projection.galerkin( -dof(uP), tf(uP), domains.omega, gauss);
+  projection.galerkin( *solution, tf(uP), domains.omega, gauss);
+  
+  msg::info << "Compute best approximation... " << msg::endl;
+  projection.pre();
+  projection.assemble();
+  projection.solve();
+
+  // save results and compute errors
+  save(+u, domains.omega, "u");
+  save(+uP, domains.omega, "uP");
+  if (abcName == "pml"){
+    save(+u, domains.omegaPml, "u_pml");
+  }
+  save(*solution, domains.omega, "u_exact");
+  save(*solution - u, domains.omega, "error");
+
+  std::complex< double > num = integrate(pow(abs(*solution - u), 2), domains.omega, gauss);
+  std::complex< double > numP = integrate(pow(abs(*solution - uP), 2), domains.omega, gauss);
+  std::complex< double > den = integrate(pow(abs(*solution), 2), domains.omega, gauss);
+  msg::info << "L_2 error = " << 100.*sqrt(num / den) << " %" << msg::endl;
+  msg::info << "Best L_2 error = " << 100.*sqrt(numP / den) << " %" << msg::endl;
+  msg::warning << "L_2 error should exclude neighbourhood elements to the point source" << msg::endl;
+  
+  for(unsigned int i = 0; i < fieldBucket.size(); ++i) {
+    delete fieldBucket[i];
+  }
+
+  return 0;
+}
diff --git a/examples/helmholtzFlow/helmholtzflow/FreeField/runTests_circle.sh b/examples/helmholtzFlow/freefield/runTests_circle.sh
similarity index 72%
rename from examples/helmholtzFlow/helmholtzflow/FreeField/runTests_circle.sh
rename to examples/helmholtzFlow/freefield/runTests_circle.sh
index 071aef80..38df525d 100644
--- a/examples/helmholtzFlow/helmholtzflow/FreeField/runTests_circle.sh
+++ b/examples/helmholtzFlow/freefield/runTests_circle.sh
@@ -6,18 +6,18 @@ do
   echo "****************************"
   echo " Test zeroth order ABC "
   echo "****************************"
-  ./demo -abcName ABC-0 -M $i -geometry 0 
+  ./example -abcName ABC-0 -M $i -geometry 0 
   echo "****************************"
   echo " Test second order ABC "
   echo "****************************"
-  ./demo -abcName ABC-2 -M $i -geometry 0
+  ./example -abcName ABC-2 -M $i -geometry 0
   echo "****************************"
   echo " Test Pade type ABC "
   echo "****************************"
-  ./demo -abcName Pade -M $i -geometry 0
+  ./example -abcName Pade -M $i -geometry 0
   echo "****************************"
   echo " Test 4 layers pml "
   echo "****************************"
-  ./demo -abcName pml -Npml 4 -M $i -geometry 0
+  ./example -abcName pml -Npml 4 -M $i -geometry 0
 done
 
diff --git a/examples/helmholtzFlow/helmholtzflow/FreeField/runTests_square.sh b/examples/helmholtzFlow/freefield/runTests_square.sh
similarity index 63%
rename from examples/helmholtzFlow/helmholtzflow/FreeField/runTests_square.sh
rename to examples/helmholtzFlow/freefield/runTests_square.sh
index 4f9d23ae..8994a558 100644
--- a/examples/helmholtzFlow/helmholtzflow/FreeField/runTests_square.sh
+++ b/examples/helmholtzFlow/freefield/runTests_square.sh
@@ -6,19 +6,19 @@ do
   echo "****************************"
   echo " Test zeroth order ABC "
   echo "****************************"
-  ./demo -abcName ABC-0 -M $i -geometry 1 
+  ./example -abcName ABC-0 -M $i -geometry 1 
   echo "****************************"
   echo " Test second order ABC "
   echo "****************************"
-  ./demo -abcName ABC-2 -M $i -geometry 1
+  ./example -abcName ABC-2 -M $i -geometry 1
   echo "****************************"
   echo " Test Pade type ABC "
   echo "****************************"
-  ./demo -abcName Pade -M $i -geometry 1 -CornerStrategy Sommerfeld
-  ./demo -abcName Pade -M $i -geometry 1 -CornerStrategy HABC
+  ./example -abcName Pade -M $i -geometry 1 -CornerStrategy Sommerfeld
+  ./example -abcName Pade -M $i -geometry 1 -CornerStrategy HABC
   echo "****************************"
   echo " Test 4 layers pml "
   echo "****************************"
-  ./demo -abcName pml -Npml 4 -M $i -geometry 1
+  ./example -abcName pml -Npml 4 -M $i -geometry 1
 done
 
diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroX/CMakeLists.txt b/examples/helmholtzFlow/helmholtzDuctHeteroX/CMakeLists.txt
deleted file mode 100644
index 4d2e4b0c..00000000
--- a/examples/helmholtzFlow/helmholtzDuctHeteroX/CMakeLists.txt
+++ /dev/null
@@ -1,8 +0,0 @@
-cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
-
-project(demo CXX)
-
-include(${CMAKE_CURRENT_SOURCE_DIR}/../../demos.cmake)
-
-add_executable(demo main.cpp ${EXTRA_INCS})
-target_link_libraries(demo ${EXTRA_LIBS})
diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroY/CMakeLists.txt b/examples/helmholtzFlow/helmholtzDuctHeteroY/CMakeLists.txt
deleted file mode 100644
index 4d2e4b0c..00000000
--- a/examples/helmholtzFlow/helmholtzDuctHeteroY/CMakeLists.txt
+++ /dev/null
@@ -1,8 +0,0 @@
-cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
-
-project(demo CXX)
-
-include(${CMAKE_CURRENT_SOURCE_DIR}/../../demos.cmake)
-
-add_executable(demo main.cpp ${EXTRA_INCS})
-target_link_libraries(demo ${EXTRA_LIBS})
diff --git a/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/CMakeLists.txt b/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/CMakeLists.txt
deleted file mode 100644
index f01d5714..00000000
--- a/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/CMakeLists.txt
+++ /dev/null
@@ -1,8 +0,0 @@
-cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
-
-project(demo CXX)
-
-include(${CMAKE_CURRENT_SOURCE_DIR}/../../../demos.cmake)
-
-add_executable(demo main.cpp ${EXTRA_INCS})
-target_link_libraries(demo ${EXTRA_LIBS})
diff --git a/examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/CMakeLists.txt b/examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/CMakeLists.txt
deleted file mode 100644
index f01d5714..00000000
--- a/examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/CMakeLists.txt
+++ /dev/null
@@ -1,8 +0,0 @@
-cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
-
-project(demo CXX)
-
-include(${CMAKE_CURRENT_SOURCE_DIR}/../../../demos.cmake)
-
-add_executable(demo main.cpp ${EXTRA_INCS})
-target_link_libraries(demo ${EXTRA_LIBS})
diff --git a/examples/helmholtzFlow/waveguide/Readme.md b/examples/helmholtzFlow/waveguide/Readme.md
new file mode 100644
index 00000000..441fbd98
--- /dev/null
+++ b/examples/helmholtzFlow/waveguide/Readme.md
@@ -0,0 +1,78 @@
+# Assessment of Local Absorbing Boundary Conditions for Heterogeneous Convected Helmholtz Problems
+> Marchner P.
+
+Dependency(ies):
+* GmshFEM : v 1.0.0
+* Gmsh : v 4.8.2
+
+## Problem description
+
+This code solves a convected Helmholtz problem in a 2d straight waveguide 
+$$\frac{D_0}{Dt}\left(\frac{1}{c_0^2}\frac{D_0}{Dt}\right) u - \rho_0^{-1}\nabla \cdot (\rho_0 \nabla ) u = 0, \quad \text{in} \; \Omega=[0, L]\times[0, H]$$
+with either a uniform (`uniform/`) or non-uniform flow (`nonUniform/`)
+
+Different ABCs are tested on the output boundary.
+For the non-uniform flow, the relative $$L^2$$-error is computed thanks to a PML reference solution.
+
+## Installation
+
+```
+  mkdir build && cd build
+  cmake ..
+  make
+```
+
+## Usage
+
+Simply run:
+```
+  ./example [PARAM]
+```
+with `[PARAM]`:
+* `-h [x]` where `[x]` is the mesh size.
+* `-L [x]` where `[x]` is the waveguide length and ABC position.
+* `-H [x]` where `[x]` is the waveguide height.
+* `-FEMorder [x]` where `[x]` is the order of the finite element hierarchical basis.
+* `-omega [x]` where `[x]` is the angular frequency for (`waveguideUniform/`) or `-kinf [x]` where `[x]` is the reference wavenumber for (`waveguideNonuniform/`).
+* `-mode [x]` where `[x]` is the mode number for the input boundary condition.
+* `-abcName [x]` where `[x]` is the ABC imposed on the right boundary.
+* `-padeOrder [x]` is the Pade order if `abcName = Pade`. 
+* `-angle [x]` is the rotation angle of the branch-cut if `abcName = Pade`.
+* `-runcase [x]` selects a specific configuration for the results reproductibility.
+
+Each subfolder contains
+* the main script (`main.cpp`)
+* a loop for testing the ABCs over a given frequency range (`main_omega.cpp`)
+To run one of the loops just change the name in the `CMakeLists.txt` and recompile
+
+## Results reproductibility
+##### Folder `uniform/`
+
+Figure 6 can be obtained with `main_omega.cpp` and
+```
+  ./example -runcase [i]
+```
+with i=1 or i=2 for respectively Figures 6a and 6b.
+
+##### Folder `nonUniform/`
+Pressure maps from Figure 8 can be obtained by compiling `main.cpp` and running
+```
+  ./example -kinf 40
+```
+The solution `vpml.msh` can be open in Gmsh.
+  
+For Figure 9, compile `main_omega.cpp` and run the command
+```
+./example -runcase 1 -L [i]
+```
+with i={0.75,1,1.25,2} to obtain the 4 sub-figures.
+
+For Figure 10, compile `main_omega.cpp` and run the command
+```
+./example -runcase 2 -L [i]
+```
+with i={0.75,1,1.25,2} to obtain the 4 sub-figures.
+
+## Reference
+
+// Preprint
diff --git a/examples/helmholtzFlow/waveguide/nonUniform/CMakeLists.txt b/examples/helmholtzFlow/waveguide/nonUniform/CMakeLists.txt
new file mode 100644
index 00000000..6f8caaab
--- /dev/null
+++ b/examples/helmholtzFlow/waveguide/nonUniform/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(example CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake)
+
+add_executable(example main.cpp ${EXTRA_INCS})
+target_link_libraries(example ${EXTRA_LIBS})
diff --git a/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/main.cpp b/examples/helmholtzFlow/waveguide/nonUniform/main.cpp
similarity index 99%
rename from examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/main.cpp
rename to examples/helmholtzFlow/waveguide/nonUniform/main.cpp
index b153a8e9..8543cdf1 100644
--- a/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/main.cpp
+++ b/examples/helmholtzFlow/waveguide/nonUniform/main.cpp
@@ -113,7 +113,6 @@ int main(int argc, char **argv)
   int FEMorder = 8;
   gmshFem.userDefinedParameter(FEMorder, "FEMorder");
   std::string gauss = "Gauss"+std::to_string(2*FEMorder+1);
-  gmshFem.userDefinedParameter(gauss, "gauss");
 
   // meshing
   double L = 2.; // duct length
@@ -279,7 +278,7 @@ int main(int argc, char **argv)
   // Compute reference solution with large PML
   msg::info << "Use a PML with N = " << Npml << " layers" << msg::endl;
   const double Lpml = Lext+Npml*h;
-  ScalarFunction< std::complex< double > > Sigma0 = 1;//4*betax*betax;
+  ScalarFunction< std::complex< double > > Sigma0 = 4;
   // Bermudez function
   ScalarFunction< std::complex< double > > SigmaX = Sigma0 / (Lpml - x< std::complex< double > >() );
   ScalarFunction< std::complex< double > > gammaX = 1-(im/k0)*SigmaX;
diff --git a/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/main_loop_omega.cpp b/examples/helmholtzFlow/waveguide/nonUniform/main_omega.cpp
similarity index 100%
rename from examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/main_loop_omega.cpp
rename to examples/helmholtzFlow/waveguide/nonUniform/main_omega.cpp
diff --git a/examples/helmholtzFlow/waveguide/uniform/CMakeLists.txt b/examples/helmholtzFlow/waveguide/uniform/CMakeLists.txt
new file mode 100644
index 00000000..6f8caaab
--- /dev/null
+++ b/examples/helmholtzFlow/waveguide/uniform/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(example CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake)
+
+add_executable(example main.cpp ${EXTRA_INCS})
+target_link_libraries(example ${EXTRA_LIBS})
diff --git a/examples/helmholtzFlow/waveguide/uniform/main.cpp b/examples/helmholtzFlow/waveguide/uniform/main.cpp
new file mode 100644
index 00000000..7522d00b
--- /dev/null
+++ b/examples/helmholtzFlow/waveguide/uniform/main.cpp
@@ -0,0 +1,266 @@
+#include "GmshFem.h"
+#include "Formulation.h"
+#include "gmsh.h"
+#include "AnalyticalFunction.h"
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::problem;
+using namespace gmshfem::domain;
+using namespace gmshfem::field;
+using namespace gmshfem::function;
+using namespace gmshfem::analytics;
+using namespace gmshfem::post;
+using namespace gmshfem::equation;
+
+void meshDuct(const double L, const double H, const double h, const int MeshOrder)
+{
+  msg::info << " - Duct of length L = " << L << msg::endl;
+  msg::info << " - Duct of heigth H = " << H << msg::endl;
+  gmsh::model::add("geometry");
+
+  // duct definition
+  gmsh::model::geo::addPoint(L, H, 0., h, 6);
+  gmsh::model::geo::addPoint(0, H, 0., h, 7);
+  gmsh::model::geo::addPoint(0, 0, 0., h, 8);
+  gmsh::model::geo::addPoint(L, 0, 0., h, 9);
+
+  gmsh::model::geo::addLine(6, 7, 5);
+  gmsh::model::geo::addLine(7, 8, 6);
+  gmsh::model::geo::addLine(8, 9, 7);
+  gmsh::model::geo::addLine(9, 6, 8);
+  gmsh::model::geo::addCurveLoop({5,6,7,8}, 1);
+  gmsh::model::geo::addPlaneSurface({1}, 1);
+  gmsh::model::geo::mesh::setTransfiniteSurface(1);
+  gmsh::model::geo::mesh::setRecombine(1,1);
+
+  // physical groups
+  gmsh::model::addPhysicalGroup(2, {1}, 1);
+  gmsh::model::setPhysicalName(2, 3, "omega");
+  gmsh::model::addPhysicalGroup(1, {5}, 2);
+  gmsh::model::setPhysicalName(1, 2, "gammaTop");
+  gmsh::model::addPhysicalGroup(1, {6}, 3);
+  gmsh::model::setPhysicalName(1, 3, "gammaLeft");
+  gmsh::model::addPhysicalGroup(1, {7}, 4);
+  gmsh::model::setPhysicalName(1, 4, "gammaBottom");
+  gmsh::model::addPhysicalGroup(1, {8}, 5);
+  gmsh::model::setPhysicalName(1, 5, "gammaRight");
+  
+  gmsh::model::addPhysicalGroup(0, {6}, 6);
+  gmsh::model::setPhysicalName(0, 6, "gammaTopCorner");
+  gmsh::model::addPhysicalGroup(0, {9}, 7);
+  gmsh::model::setPhysicalName(0, 7, "gammaBottomCorner");
+
+  // generate mesh
+  gmsh::model::geo::synchronize();
+  gmsh::model::mesh::generate(2);
+  gmsh::model::mesh::setOrder(MeshOrder);
+  gmsh::option::setNumber("Mesh.ElementOrder",MeshOrder);
+  gmsh::model::mesh::setOrder(MeshOrder);
+}
+
+int main(int argc, char **argv)
+{
+  //*****
+  // Problem declaration
+  //*****
+  GmshFem gmshFem(argc, argv);
+  const double pi = 3.14159265358979323846264338327950288;
+  const std::complex< double > im(0., 1.);
+  
+  // numerical parameters
+  int FEMorder = 8;
+  gmshFem.userDefinedParameter(FEMorder, "FEMorder");
+  std::string gauss = "Gauss"+std::to_string(2*FEMorder+1);
+
+  // meshing
+  double L = 0.5; // duct length
+  gmshFem.userDefinedParameter(L, "L");
+  double H = 0.25; // duct height
+  gmshFem.userDefinedParameter(H, "H");
+  double h = 1./40; // meshsize
+  gmshFem.userDefinedParameter(h, "h");
+  meshDuct(L, H, h, 1);
+ 
+  // duct problem specifications
+  std::string problem = "hard";
+  gmshFem.userDefinedParameter(problem, "problem");
+  double w = 30; // frequency
+  gmshFem.userDefinedParameter(w, "omega");
+  std::vector< unsigned int > n = {3}; // single mode
+  // std::vector< unsigned int > n = {0,2,4,6,8}; // multiple modes modes
+  std::vector< double > A(n.size(),1.); // select amplitudes
+  int N_modes = n.size();
+  msg::info << " There are " << N_modes << " input modes " << msg::endl;
+ 
+  // Choose ABC type
+  std::string abcName = "Pade";
+  gmshFem.userDefinedParameter(abcName, "abcName");
+
+  // mean flow
+  const double M = 0.8;
+  double beta = sqrt(1-M*M);
+
+  std::vector< std::complex< double > > KX(N_modes);
+  std::vector< double > KY(N_modes), SIGMA(N_modes);
+  for(unsigned int j = 0; j < n.size(); ++j) {
+    KY[j] = n[j]*pi/H;
+    SIGMA[j] = w*w-beta*beta*KY[j]*KY[j];
+    if ( SIGMA[j] >=0 ) {
+      KX[j] = (1/(beta*beta))*(-M*w + sqrt(SIGMA[j]));
+      msg::info << " Mode " << n[j] << " is propagative " << msg::endl;
+      if ( w>(beta*KY[j]) && w<KY[j] ) {
+        msg::info << " Mode " << n[j] << " is inverse upstream ! " << msg::endl;
+      }
+    }
+    else {
+      KX[j] = (1/(beta*beta))*(-M*w - im* sqrt(abs(SIGMA[j])));
+      msg::info << " Mode " << n[j] << " is evanescent " << msg::endl;
+    }
+    msg::info << " - propagating wavenumber, kx = " << KX[j] << " of mode " << n[j] << msg::endl;
+  }
+
+  float PtsByWl = 2*pi*FEMorder*(1+M) / w / h;
+  msg::info << "********************************" << msg::endl;
+  msg::info << " - convected Helmholtz duct problem " << msg::endl;
+  msg::info << " - wave number = " << w << msg::endl;
+  msg::info << " - FEM basis order = " << FEMorder << "" << msg::endl;
+  msg::info << " - Approximate dofs by wavelength = " << PtsByWl << "" << msg::endl;
+  msg::info << "********************************" << msg::endl;
+
+  // declare formulation and FEM domains
+  std::vector< FieldInterface< std::complex< double > > * > fieldBucket;
+  Formulation< std::complex< double > > formulation("helmholtzflow");
+
+  Domain omega(2, 1);
+  Domain gammaTop(1, 2);
+  Domain gammaLeft(1, 3);
+  Domain gammaBottom(1, 4);
+  Domain gammaRight(1, 5);
+  Domain gammaTopCorner(0, 6); 
+  Domain gammaBottomCorner(0, 7); 
+
+  Field< std::complex< double >, Form::Form0 > v("v", omega | gammaLeft | gammaRight | gammaTop | gammaBottom, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder);
+  // define analytical solution
+  Function< std::complex< double >, Degree::Degree0 > *solution = nullptr;
+
+  if(problem == "soft") {
+    v.addConstraint(gammaTop, 0.);
+    v.addConstraint(gammaBottom, 0.);
+    solution = new AnalyticalFunction< helmholtz2D::DuctModeSolutionMulti< std::complex< double > > >(w, M, H, A, n, 0., 0., 1);
+  }
+  else {
+    solution = new AnalyticalFunction< helmholtz2D::DuctModeSolutionMulti< std::complex< double > > >(w, M, H, A, n, 0., 0., 0);
+  }
+  solution->activateMemorization();
+
+  formulation.galerkin(vector< std::complex< double > >(beta*beta, 0., 0.) * grad(dof(v)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omega, gauss);
+  formulation.galerkin(vector< std::complex< double > >(0., 1., 0.) * grad(dof(v)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(v)), omega, gauss);
+  formulation.galerkin(vector< std::complex< double > >(im*w*M, 0., 0.) * grad(dof(v)), tf(v), omega, gauss);
+  formulation.galerkin(-im*w*M*dof(v), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omega, gauss);
+  formulation.galerkin(- w * w * dof(v), tf(v), omega, gauss);
+
+  // input BC
+  ScalarFunction< std::complex< double > > Input_modal_sum=0.;
+  for(unsigned int j = 0; j < n.size(); ++j) {
+    if(problem == "soft"){
+      Input_modal_sum = Input_modal_sum + A[j] * KX[j] * sin( KY[j] * y< std::complex< double > >() );
+    }
+    else if(problem == "hard"){
+      Input_modal_sum = Input_modal_sum + A[j] * KX[j] * cos( KY[j] * y< std::complex< double > >() );
+    }
+    else {
+      msg::error << " this problem is not defined " << msg::endl;
+    }
+  }
+  formulation.galerkin(-beta*beta * im * Input_modal_sum , tf(v), gammaLeft, gauss);
+  formulation.galerkin(-im * w * M * dof(v), tf(v), gammaLeft, gauss);
+
+  // output BC
+  if(abcName == "ABC-0") {
+    msg::info << "Use Sommerfeld ABC." << msg::endl;
+    const double angle = 0.;
+    std::complex< double > wM = w*(M-exp(im*angle/2.)) / (beta*beta); //angle=0 -> -k/(1+M);
+    formulation.galerkin(-beta*beta * im * wM * dof(v), tf(v), gammaRight, gauss);
+    formulation.galerkin(im * w * M * dof(v), tf(v), gammaRight, gauss);
+  }
+  else if(abcName == "ABC-2") {
+    msg::info << "Use second order ABC." << msg::endl;
+    const double angle = 0.;
+    // contribution to the mass and rigidity matrices
+    std::complex< double > wM = w*(M-cos(angle/2.)) / (beta*beta);
+    std::complex< double > wK = -exp(-im*angle/2.)/(2*w);
+    formulation.galerkin(-beta*beta * im * wM * dof(v), tf(v), gammaRight, gauss);
+    formulation.galerkin( beta*beta * im * wK * grad(dof(v)), grad(tf(v)), gammaRight, gauss);
+    formulation.galerkin(im * w * M * dof(v), tf(v), gammaRight, gauss);
+    
+    /* use DtN eigenvalues instead of Laplace Beltrami operator
+    std::complex< double > kT = -k + KY[0]*KY[0]/(2*k);
+    formulation.galerkin(-beta*beta * im * kT * dof(v), tf(v), gammaRight, gauss);
+    */
+  }
+  else if(abcName == "Pade") {
+    int padeOrder = 4;
+    gmshFem.userDefinedParameter(padeOrder, "padeOrder");
+    double angle = -pi/4.;
+    gmshFem.userDefinedParameter(angle, "angle");
+    msg::info << "Use Pade ABC of order " << padeOrder << " with angle " << angle << " rad"<< msg::endl;
+
+    const double Np = 2. * padeOrder + 1.;
+    const std::complex< double > exp1 = std::complex<double>(std::cos(angle),std::sin(angle));
+    const std::complex< double > exp2 = std::complex<double>(std::cos(angle/2.),std::sin(angle/2.));
+    const std::complex< double > coef = 2./Np;
+    std::vector< std::complex< double > > c(padeOrder, 0.);
+    for(int i = 0; i < padeOrder; ++i) {
+      c[i] = std::tan((i + 1) * pi / Np);
+      c[i] *= c[i];
+    }
+
+    // define the auxiliary fields
+    std::vector< Field< std::complex< double >, Form::Form0 >* > phi;
+    Field< std::complex< double >, Form::Form0 >* phis; 
+    for(int i = 0; i < padeOrder; ++i) {
+      phis = new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(i), gammaRight | gammaTopCorner | gammaBottomCorner, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder);
+      if(problem == "soft") { // impose Dirichlet constraints on the auxiliary fields
+        phis->addConstraint(gammaTopCorner,0.);
+        phis->addConstraint(gammaBottomCorner,0.);
+      }
+      phi.push_back(phis);
+      fieldBucket.push_back(phi.back());
+    }
+    
+    // write the augmented weak form
+    formulation.galerkin(im * w * exp2 * dof(v), tf(v), gammaRight, gauss);
+    for(int i = 0; i < padeOrder; ++i) {
+      // boundary integral terms relating the auxiliary fields
+      formulation.galerkin(im * w * exp2 * coef * c[i] * dof(*phi[i]), tf(v), gammaRight, gauss);
+      formulation.galerkin(im * w * exp2 * coef * c[i] * dof(v), tf(v), gammaRight, gauss);
+
+      // coupling of the auxiliary equations
+      formulation.galerkin( grad(dof(*phi[i])), grad(tf(*phi[i])), gammaRight, gauss);
+      formulation.galerkin(- (w*w)/(beta*beta) * (exp1 * c[i] + 1.) * dof(*phi[i]), tf(*phi[i]), gammaRight, gauss);
+      formulation.galerkin(- (w*w)/(beta*beta) * exp1 * (c[i] + 1.) * dof(v), tf(*phi[i]), gammaRight, gauss);
+    }
+  }
+  else if(abcName == "DtN") {
+    msg::info << "Use analytical DtN." << msg::endl;
+    formulation.galerkin(beta*beta * im * KX[0] * dof(v), tf(v), gammaRight, gauss);
+    formulation.galerkin(im * w * M * dof(v), tf(v), gammaRight, gauss);
+  }
+
+  formulation.pre();
+  formulation.assemble();
+  formulation.solve();
+  save(+v, omega, "v");
+  save(*solution, omega, "v_exact");
+  save(*solution - v, omega, "error");
+  std::complex< double > num = integrate(pow(abs(*solution - v), 2), omega, gauss);
+  std::complex< double > den = integrate(pow(abs(*solution), 2), omega, gauss);
+  msg::info << "L_2 error = " << 100.*sqrt(num / den) << " %" << msg::endl;
+
+  for(unsigned int i = 0; i < fieldBucket.size(); ++i) {
+    delete fieldBucket[i];
+  }
+
+  return 0;
+}
diff --git a/examples/helmholtzFlow/helmholtzflow/Duct_ABCs/main_omega.cpp b/examples/helmholtzFlow/waveguide/uniform/main_omega.cpp
similarity index 100%
rename from examples/helmholtzFlow/helmholtzflow/Duct_ABCs/main_omega.cpp
rename to examples/helmholtzFlow/waveguide/uniform/main_omega.cpp
-- 
GitLab