diff --git a/examples/helmholtzFlow/freefield/README.md b/examples/helmholtzFlow/freefield/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..10fa70050882165ad47d7a7d7f2ac6189c0bb66d
--- /dev/null
+++ b/examples/helmholtzFlow/freefield/README.md
@@ -0,0 +1,39 @@
+# Absorbing boundary condtions for free field convected propagation
+
+Dependency(ies):
+* GmshFEM : v 1.0.0
+* Gmsh : v 4.8.2
+
+## Problem description
+
+These routines solve a convected Helmholtz propagation problem with a Dirac right-hand-side, either in a circular or square domain. 
+Different ABCs are tested on the output boundary, and are compared to the analytical solution.
+
+## Installation and usage 
+Simply run 
+
+```
+  mkdir build && cd build
+  cmake ..
+  make
+  ./example [PARAM]
+```
+with `[PARAM]`:
+* `-geometry [x]` where `[x]` is the domain geometry (0: disk, 1: square)
+* `-k [x]` is the freefield wavenumber, $`k>0`$
+* `-R [x]` where `[x]` is the disk radius or square size, $`R>0`$
+* `-xs [x]` and `-ys [x]` specify the point source location
+* `-M [x]` is the mean flow Mach number, $`0<M<1`$
+* `-theta [x]` is the mean flow angle orientation, $\theta \in [0, 2\pi]$
+* `-abcName [x]` is the ABC type. The available choices are `ABC-0`, `ABC-2`, `Pade` and `pml`
+* `-padeOrder [x]` and `-angle [x]` specify respectively the number of auxiliary function and branch rotation angle for the `Pade` abc.
+
+By default a circular domain with a Padé condition is used, for a centered source at the frequency $k=6\pi$, with a mean flow at $M=0.8$ oriented at $\theta=\pi/4$.
+
+## Reference
+> To be updated
+
+## Results reproducibility
+We assess the performance of different ABCs in the case of a radiating point source set in a circular domain.
+A loop on the Mach number is launched by executing the file `runTests_ABC_circle.sh`, which allows to obtain the data shown at the end of Section 2.
+
diff --git a/examples/helmholtzFlow/freefield/main.cpp b/examples/helmholtzFlow/freefield/main.cpp
index c6d5bfe5ca5c345057ff2162197c7159a6e9084d..2f5e3b56a5eaefcc1761cdfaeba28785ab926f28 100644
--- a/examples/helmholtzFlow/freefield/main.cpp
+++ b/examples/helmholtzFlow/freefield/main.cpp
@@ -246,10 +246,10 @@ int main(int argc, char **argv)
   // geometry parameters
   int geometry = 0; // 0-circle, 1-square
   gmshFem.userDefinedParameter(geometry, "geometry");
-  double R = 0.5; // Circle radius or square of size [-R,R] x [-R,R]
+  double R = 1.; // Circle radius or square of size [-R,R] x [-R,R]
   gmshFem.userDefinedParameter(R, "R");
   // Mesh order
-  int ElemOrder = 6;
+  int ElemOrder = 4;
   // Point source location
   double xs = 0.;
   double ys = 0.;
@@ -259,26 +259,40 @@ int main(int argc, char **argv)
 
   // physical parameters
   const double pi = 3.14159265358979323846264338327950288;
-  const double k = 6 * pi; // free field wavenumber
+  double k = 6 * pi; // free field wavenumber
+  gmshFem.userDefinedParameter(k, "k");
   const std::complex< double > im(0., 1.);
   double M = 0.8; // 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 = 5; // FEM shape function order
-  std::string gauss = "Gauss" + std::to_string(2 * order + 1); // increase Nr of Gauss points when order > 6
-  double lc = 0.02; // mesh size
-  gmshFem.userDefinedParameter(lc, "lc");
-
+  int order; // FEM shape function order
+  if (M >= 0.9) {
+    order = 9;
+  }
+  else if (M < 0.9 && M >= 0.7) {
+    order = 6;
+  }
+  else {
+    order = 4;
+  }
+  
+  std::string gauss = "Gauss" + std::to_string(2 * order + 1);
+  double pointsByWl = 8; // pointsByWl = (2 * pi * order * (1 - abs(M))) / (lc * k);
+  double lc = (2 * pi * order * (1 - abs(M))) / (pointsByWl * k); // choose lc based on pointsByWl
+  msg::info << " - pointsByWl =  " << (2 * pi * order * (1 - abs(M))) / (lc * k) << msg::endl;
+  msg::info << " - lc =  " << lc << msg::endl;
+  msg::info << " - FEM basis order = " << order << "" << msg::endl;
+ 
   // Choose ABC - exterior boundary condition
   std::string abcName = "Pade"; // available choices: ABC-0, ABC-2, Pade, pml
   gmshFem.userDefinedParameter(abcName, "abcName");
-  bool L0 = 1; // zeroth-order symbol
+  bool L0 = false; // activate zeroth-order symbol
   gmshFem.userDefinedParameter(L0, "L0");
-  bool withCorner = true;
+  bool withCorner = true; // try corner treatment for square domain 
   gmshFem.userDefinedParameter(withCorner, "withCorner");
   int Npml;
   double Rpml;
@@ -303,13 +317,6 @@ int main(int argc, char **argv)
     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);
@@ -354,6 +361,7 @@ int main(int argc, char **argv)
     formulation.integral(im * k * (1 - Mn) * dof(u), tf(u), domains.gammaPhy, gauss);
     formulation.integral(-Mn * Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss);
     if ( L0 && (geometry==0) ) {
+      msg::info << "Use curvature correction" << msg::endl;
       formulation.integral( + (beta * 0.5 / R) * dof(u), tf(u), domains.gammaPhy, gauss);
     }
   }
@@ -374,7 +382,7 @@ int main(int argc, char **argv)
     }
     if ( withCorner && (geometry==1) ) {
       msg::info << "with corner treatment" << msg::endl;
-      // corner condition 2nd order abc, M = 0
+      // try a corner condition for 2nd order abc, M = 0
       for(int ci = 0; ci < 4; ++ci) {
         formulation.integral( 0.75 * dof(u), tf(u), domains.corner[ci], gauss);
       }
@@ -419,7 +427,7 @@ int main(int argc, char **argv)
         formulation.integral(- (beta * beta) * grad(dof(*phi[i])), grad(tf(*phi[i])), domains.gammaPhy, gauss);
         formulation.integral(-2. * im * k * Mt * tangent< std::complex< double > >() * grad(dof(*phi[i])), tf(*phi[i]), domains.gammaPhy, gauss);
         // regularization
-        std::complex<double> keps = k-0.*im*0.4*pow(k,(1/3.))*pow((beta/R),2/3.);
+        std::complex<double> keps = k;//-0.*im*0.4*pow(k,(1/3.))*pow((beta/R),2/3.);
         formulation.integral((keps * keps) * (exp1 * c[i] + 1.) * dof(*phi[i]), tf(*phi[i]), domains.gammaPhy, gauss);
         formulation.integral((keps * keps) * exp1 * (c[i] + 1.) * dof(u), tf(*phi[i]), domains.gammaPhy, gauss);
       }
@@ -608,12 +616,12 @@ int main(int argc, char **argv)
   std::complex< double > num = integrate(pow(abs(*solution - u), 2), domains.omegaErr, gauss);
   std::complex< double > numP = integrate(pow(abs(*solution - uP), 2), domains.omegaErr, gauss);
   std::complex< double > den = integrate(pow(abs(*solution), 2), domains.omegaErr, gauss);
-  msg::info << "L_2 error = " << 100. * sqrt(num / den) << " %" << msg::endl;
+  msg::info << "Relative L_2 error = " << 100. * sqrt(num / den) << " %" << msg::endl;
   msg::info << "Best L_2 error = " << 100. * sqrt(numP / den) << " %" << msg::endl;
 
   std::complex< double > numGamma = integrate(pow(abs(*solution - u), 2), domains.gammaPhy, gauss);
   std::complex< double > denGamma = integrate(pow(abs(*solution), 2), domains.gammaPhy, gauss);
-  msg::info << "L_2 error on interface = " << 100. * sqrt(numGamma / denGamma) << " %" << msg::endl;
+  msg::info << "Interface relative L_2 error = " << 100. * sqrt(numGamma / denGamma) << " %" << msg::endl;
 
   for(unsigned int i = 0; i < fieldBucket.size(); ++i) {
     delete fieldBucket[i];
diff --git a/examples/helmholtzFlow/freefield/runTests_ABC_circle.sh b/examples/helmholtzFlow/freefield/runTests_ABC_circle.sh
new file mode 100644
index 0000000000000000000000000000000000000000..aa8495d6d90a12eb1515f0fec04b4a76b431b228
--- /dev/null
+++ b/examples/helmholtzFlow/freefield/runTests_ABC_circle.sh
@@ -0,0 +1,32 @@
+cd build
+cmake .. && make
+
+for i in 0 0.05 0.1 #0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95
+do
+  echo "---------- M = $i ----------"
+  echo " Test zeroth order ABC "
+  echo "****************************"
+  ./example -abcName ABC-0 -M $i -geometry 0 > res.txt
+  ERR=$(grep "Relative L_2 error =" res.txt | awk '{print $7}')
+  echo "Relative L2 error $ERR" 
+  echo $ERR >> ABC0.txt
+  echo " Test second order ABC "
+  echo "****************************"
+  ./example -abcName ABC-2 -M $i -geometry 0 > res.txt
+  ERR=$(grep "Relative L_2 error =" res.txt | awk '{print $7}')
+  echo "Relative L2 error $ERR" 
+  echo $ERR >> ABC2.txt
+  echo " Test Pade order 4 ABC "
+  echo "****************************"
+  ./example -abcName Pade -M $i -geometry 0 -padeOrder 4 > res.txt
+  ERR=$(grep "Relative L_2 error =" res.txt | awk '{print $7}') 
+  echo "Relative L2 error $ERR" 
+  echo $ERR >> Pade.txt
+  BEST=$(grep "Best L_2 error =" res.txt | awk '{print $7}')
+  echo "Best L2 error $BEST"
+  echo $BEST >> Best.txt
+  echo $i >> Mach.txt
+done
+
+paste Mach.txt ABC0.txt ABC2.txt Pade.txt Best.txt > ABC_results.txt
+rm Mach.txt ABC0.txt ABC2.txt Pade.txt Best.txt
\ No newline at end of file
diff --git a/examples/helmholtzFlow/hotjet/CMakeLists.txt b/examples/helmholtzFlow/hotjet/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e9abdae9076c426ebcf8a342e2e9cf81cb5be22a
--- /dev/null
+++ b/examples/helmholtzFlow/hotjet/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 mesh.cpp flowPmlFunctions.cpp ${EXTRA_INCS})
+target_link_libraries(example ${EXTRA_LIBS})
\ No newline at end of file
diff --git a/examples/helmholtzFlow/hotjet/flowPmlFunctions.cpp b/examples/helmholtzFlow/hotjet/flowPmlFunctions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..099b362add2f9b2ce4c9e0e3568f74484ac41ec9
--- /dev/null
+++ b/examples/helmholtzFlow/hotjet/flowPmlFunctions.cpp
@@ -0,0 +1,37 @@
+#include "flowPmlFunctions.h"
+
+TensorFunction< std::complex< double > > getInverseLorentzTensor(
+    ScalarFunction< std::complex< double > > Mx, ScalarFunction< std::complex< double > > My
+    , ScalarFunction< std::complex< double > > beta)
+{
+  ScalarFunction< std::complex< double > > Alphax = 1 + Mx * Mx / (beta * (1 + beta));
+  ScalarFunction< std::complex< double > > Alphay = 1 + My * My / (beta * (1 + beta));
+  ScalarFunction< std::complex< double > > K = Mx * My / (beta * (1 + beta));
+  return tensor< std::complex< double > >(beta * Alphay, -beta * K, 0., -beta * K, beta * Alphax, 0., 0., 0., 0.);
+}
+
+TensorFunction< std::complex< double > > getInversePMLJacobian(
+    ScalarFunction< std::complex< double > > k, std::pair<double, double> xlim, std::pair<double, double> ylim, double Lpml,
+    ScalarFunction< std::complex< double > > beta,
+    ScalarFunction< std::complex< double > > &detJ)
+{
+  std::complex< double > im(0., 1.);
+  ScalarFunction< std::complex< double > > Sigma0 = beta;
+  // x-PML
+  double Lxpml1 = xlim.first - Lpml; // PML left limit in x (-)
+  double Lxpml2 = xlim.second + Lpml; // PML right limit in x (+)
+  ScalarFunction< std::complex< double > > SigmaX = heaviside(-(x< std::complex< double > >() - xlim.first) ) * Sigma0 / ( abs(Lxpml1) - abs(x< std::complex< double > >()) )
+    + heaviside( x< std::complex< double > >() - xlim.second ) * Sigma0 / ( abs(Lxpml2) - abs(x< std::complex< double > >()) );
+  
+  // y-PML
+  double Lypml1 = ylim.first - Lpml; // PML lower limit in y (-)
+  double Lypml2 = ylim.second + Lpml; // PML upper limit in y (+) 
+  ScalarFunction< std::complex< double > > SigmaY = heaviside(-(y< std::complex< double > >() - ylim.first) ) * Sigma0 / ( abs(Lypml1) - abs(y< std::complex< double > >()) )
+    + heaviside( y< std::complex< double > >() - ylim.second ) * Sigma0 / ( abs(Lypml2) - abs(y< std::complex< double > >()) );
+  // PML stretchings
+  ScalarFunction< std::complex< double > > gammaX = 1 - (im / k) * SigmaX;
+  ScalarFunction< std::complex< double > > gammaY = 1 - (im / k) * SigmaY;
+
+  detJ = gammaX * gammaY;
+  return tensor< std::complex< double > >(1. / gammaX, 0., 0., 0., 1. / gammaY, 0., 0., 0., 0.);
+}
\ No newline at end of file
diff --git a/examples/helmholtzFlow/hotjet/flowPmlFunctions.h b/examples/helmholtzFlow/hotjet/flowPmlFunctions.h
new file mode 100644
index 0000000000000000000000000000000000000000..2e8999adcd255dc335112edce3985f544f3aa58f
--- /dev/null
+++ b/examples/helmholtzFlow/hotjet/flowPmlFunctions.h
@@ -0,0 +1,16 @@
+#include <gmshfem/Formulation.h>
+#include <gmshfem/GmshFem.h>
+
+using namespace gmshfem;
+using namespace gmshfem::function;
+
+TensorFunction< std::complex< double > > getInverseLorentzTensor(
+    ScalarFunction< std::complex< double > > Mx, 
+    ScalarFunction< std::complex< double > > My, 
+    ScalarFunction< std::complex< double > > beta);
+
+TensorFunction< std::complex< double > > getInversePMLJacobian(
+    ScalarFunction< std::complex< double > > k, 
+    std::pair<double, double> xlim, std::pair<double, double> ylim, double Lpml,
+    ScalarFunction< std::complex< double > > beta,
+    ScalarFunction< std::complex< double > > &detJ);
\ No newline at end of file
diff --git a/examples/helmholtzFlow/hotjet/main.cpp b/examples/helmholtzFlow/hotjet/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d7938ddcc90ee101a61bd0a62bda8591ed8fbbdd
--- /dev/null
+++ b/examples/helmholtzFlow/hotjet/main.cpp
@@ -0,0 +1,142 @@
+//#include <gmshfem/Formulation.h>
+//#include <gmshfem/GmshFem.h>
+#include "mesh.h"
+#include "flowPmlFunctions.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::post;
+using namespace gmshfem::equation;
+
+int main(int argc, char **argv)
+{
+  GmshFem gmshFem(argc, argv);
+  
+  // 
+  double pi = 3.14159265358979323846264338327950288;
+  std::complex< double > im(0., 1.);
+  std::string abcType = "pml"; // available choices: sommerfeld, pml
+  gmshFem.userDefinedParameter(abcType, "abcType");
+
+  // 
+  double Sigma = 1.0018;
+  gmshFem.userDefinedParameter(Sigma, "Sigma");
+  int ElemOrder = 2;
+  int Npml = 4;
+  gmshFem.userDefinedParameter(Npml, "Npml");
+  int RefSource = 18;
+
+  std::pair<double, double> xlim;
+  std::pair<double, double> ylim;
+  double lc;
+  meshJet(Sigma, RefSource, Npml, ElemOrder, xlim, ylim, lc);
+  
+  msg::info << "mesh size: " << lc << msg::endl;
+  Domain omegaPhy = Domain("omega");
+  Domain gammaABC = Domain("gamma");
+  Domain omegaPml = Domain("omega_pml");
+  Domain omega;
+  if (abcType == "pml") {
+    omega = omegaPhy | omegaPml;
+    msg::info << "x-pml-limit-left: " << xlim.first << msg::endl;
+    msg::info << "x-pml-limit-right: " << xlim.second << msg::endl;
+    msg::info << "y-pml-limit-down: " << ylim.first << msg::endl;
+    msg::info << "y-pml-limit-up: " << ylim.second << msg::endl;
+  }
+  else {
+    omega = omegaPhy;
+  }
+
+  // flow acoustic problem
+  double w = 200*pi;
+  int FEMorder = 5;
+  std::string gauss = "Gauss" + std::to_string(2 * FEMorder + 1);
+  double SigmaS = 0.136*Sigma;
+  
+  // define source
+  ScalarFunction< std::complex< double > > source = SigmaS*exp(- r2d<std::complex< double > >() / (2.*SigmaS*SigmaS) );
+  // save(source, omega, "source");
+  // define mean flow
+  double Tj = 600; // Jet temperature [Kelvin]
+  double gamma = 1.4;
+  double R = 287.;
+  double Mj = 0.9;
+  double cj = sqrt(gamma*R*Tj);
+  double uj = Mj*cj;
+  ScalarFunction< std::complex< double > > u0y = uj*exp(- y<std::complex< double > >()*y<std::complex< double > >() / (2.*Sigma*Sigma) );
+  
+  double Tinf = 300; // Reference Temperature [Kelvin]
+  double p0 = 103330;
+  double rhoj = gamma*p0/(cj*cj);
+  ScalarFunction< std::complex< double > > rho0_inv = Tinf/Tj - (Tinf/Tj - 1.) * (u0y/uj) + ((gamma-1)/2) * (Mj*Mj) * (u0y/uj) * (1-u0y/uj);
+  ScalarFunction< std::complex< double > > rho0 = rhoj/rho0_inv;
+  ScalarFunction< std::complex< double > > c0 = sqrt(gamma*p0/rho0);
+   
+  ScalarFunction< std::complex< double > > Mx = u0y/c0;
+  ScalarFunction< std::complex< double > > My = 0.;
+  ScalarFunction< std::complex< double > > k = w/c0;
+
+  save(rho0, omega, "rho0");
+  save(u0y, omega, "u0");
+  save(u0y/c0, omega, "Mach");
+
+  // set variational formulation
+  Formulation< std::complex< double > > formulation("helmholtzflow");
+  Field< std::complex< double >, Form::Form0 > u("u", omega | gammaABC, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder);
+
+  // convected Helmholz weak form
+  formulation.integral((1./rho0) * vector< std::complex< double > >(1 - Mx * Mx, -Mx * My, 0.) * grad(dof(u)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(u)), omegaPhy, gauss);
+  formulation.integral((1./rho0) * vector< std::complex< double > >(-Mx * My, 1 - My * My, 0.) * grad(dof(u)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), omegaPhy, gauss);
+  formulation.integral(- (1./rho0) * (k*k) * dof(u), tf(u), omegaPhy, gauss);
+  formulation.integral((1./rho0) * dof(u), vector< std::complex< double > >(-im * k * Mx, -im * k * My, 0.) * grad(tf(u)), omegaPhy, gauss, term::ProductType::Scalar);
+  formulation.integral((1./rho0) * vector< std::complex< double > >(im * k * Mx, im * k * My, 0.) * grad(dof(u)), tf(u), omegaPhy, gauss);
+
+  formulation.integral(source, tf(u), omegaPhy, gauss);
+
+  
+  if (abcType == "sommerfeld") { // Zeroth-order ABC
+    msg::info << "use sommerfeld abc " << msg::endl;
+    formulation.integral((1./rho0) * im * k * dof(u), tf(u), gammaABC, gauss);
+  }
+  else if(abcType == "pml") { // PML
+    msg::info << "use a pml " << msg::endl;
+    ScalarFunction< std::complex< double > > beta = sqrt(1-Mx*Mx-My*My);
+    TensorFunction< std::complex< double > > Linv;
+    Linv = getInverseLorentzTensor(Mx, My, beta);
+
+    ScalarFunction< std::complex< double > > detJ;
+    TensorFunction< std::complex< double > > JPml_invT;
+    JPml_invT = getInversePMLJacobian(k, xlim, ylim, Npml*lc, beta, detJ);
+
+    VectorFunction< std::complex< double > > MM = vector< std::complex< double > >(Mx, My, 0.);
+    VectorFunction< std::complex< double > > JPml_invT_M = JPml_invT * MM;
+    TensorFunction< std::complex< double > > JPml_Linv = JPml_invT * Linv;
+
+    // stabilized PML variational formulation
+    formulation.integral((1./rho0) * detJ * JPml_Linv*grad(dof(u)) , JPml_Linv*grad(tf(u)) , omegaPml, gauss, term::ProductType::Scalar);
+    formulation.integral((1./rho0) * k * k/(beta * beta) * detJ * JPml_invT_M * dof(u) , JPml_invT_M * tf(u), omegaPml, gauss, term::ProductType::Scalar);
+    formulation.integral((1./rho0) * im * k / beta * detJ * JPml_Linv * grad(dof(u)), JPml_invT_M * tf(u), omegaPml, gauss, term::ProductType::Scalar);
+    formulation.integral(- (1./rho0) * im * k /beta * detJ * JPml_invT_M * dof(u), JPml_Linv * grad(tf(u)), omegaPml, gauss, term::ProductType::Scalar);
+    formulation.integral(- (1./rho0) * k * k/(beta * beta) * detJ * dof(u), tf(u), omegaPml, gauss);
+  }
+  else {
+    msg::error << "ABC Type not available !" << msg::endl;
+    exit(0);
+  }
+
+
+  // solve
+  formulation.pre();
+  formulation.assemble();
+  formulation.solve();
+
+  // save
+  save(+u, omega, "u");
+
+
+  return 0;
+}
\ No newline at end of file
diff --git a/examples/helmholtzFlow/hotjet/mesh.cpp b/examples/helmholtzFlow/hotjet/mesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9c2e10d9eba9f6c07e1de57b7b510ebb0bd8a9a3
--- /dev/null
+++ b/examples/helmholtzFlow/hotjet/mesh.cpp
@@ -0,0 +1,189 @@
+#include "mesh.h"
+
+void meshJet(const double sigma, const double SourceMeshRefinement, const int NrOfLayers, const int MeshElemOrder, std::pair<double, double> &xlim, std::pair<double, double> &ylim, double &lc)
+{
+  // gmsh::initialize();
+  std::cout << " - Sigma = " << sigma << std::endl;
+  std::cout << " - Number of PML layers = " << NrOfLayers << std::endl;
+  std::cout << " - Source refinement factor = " << SourceMeshRefinement << std::endl;
+  gmsh::model::add("geometry");
+  
+  // parameters
+  lc = 0.6*sigma;
+  double Lx = 7.*sigma;
+  double Ly = 6.*sigma;
+
+  double Cy = 3.*sigma;
+  double Lx_down = 9.*sigma;
+
+  int Nx = std::ceil((2.0*Lx+Lx_down)/lc);
+  std::cout << "Nx = " << Nx << std::endl;
+  // int Ny_long = std::ceil(Ly/lc);
+  int Ny_short = std::ceil(((Ly-Cy))/lc);
+  std::cout << "Ny_short = " << Ny_short << std::endl;
+  double Lpml = NrOfLayers*lc;
+  double MeshSizeSource = lc/SourceMeshRefinement;
+
+  // geometry definition
+  gmsh::model::geo::addPoint(0., 0., 0., MeshSizeSource, 1);
+  gmsh::model::geo::addPoint(-4*sigma, 0., 0., MeshSizeSource, 501);
+  gmsh::model::geo::addPoint(-3*sigma, 0., 0., MeshSizeSource, 502);
+  gmsh::model::geo::addPoint(-2*sigma, 0., 0., MeshSizeSource, 503);
+  for (int j=0; j<=15 ; j++) {
+    gmsh::model::geo::addPoint((j+2)*sigma, 0., 0., MeshSizeSource, j+504);
+  }
+  gmsh::model::geo::addPoint(-sigma, 0., 0., MeshSizeSource, 520);
+  gmsh::model::geo::addPoint(sigma, 0., 0., MeshSizeSource, 521);
+  
+  gmsh::model::geo::addPoint(-Lx, -Cy, 0., MeshSizeSource, 523);
+  gmsh::model::geo::addPoint(Lx+Lx_down, -Cy, 0., MeshSizeSource, 524);
+ 
+  gmsh::model::geo::addPoint(-Lx, -Ly, 0., lc, 2);
+  gmsh::model::geo::addPoint(Lx+Lx_down, -Ly, 0, lc, 3);
+  gmsh::model::geo::addPoint(Lx+Lx_down, Ly, 0, lc, 4);
+  gmsh::model::geo::addPoint(-Lx, Ly, 0, lc, 5);
+
+  gmsh::model::geo::addPoint(-Lx-Lpml, -Ly, 0, lc, 6);
+  gmsh::model::geo::addPoint(-Lx, -Ly-Lpml, 0, lc, 7);
+  gmsh::model::geo::addPoint(Lx+Lpml+Lx_down, -Ly, 0, lc, 8);
+  gmsh::model::geo::addPoint(Lx+Lx_down, -Ly-Lpml, 0, lc, 9);
+  
+  gmsh::model::geo::addPoint(Lx+Lpml+Lx_down, Ly, 0, lc, 10);
+  gmsh::model::geo::addPoint(Lx+Lx_down, Ly+Lpml, 0, lc, 11);
+  gmsh::model::geo::addPoint(-Lx-Lpml, Ly, 0, lc, 12);
+  gmsh::model::geo::addPoint(-Lx, Ly+Lpml, 0, lc, 13);
+
+  gmsh::model::geo::addPoint(-Lx-Lpml, -Ly-Lpml, 0, lc, 14);
+  gmsh::model::geo::addPoint(Lx+Lpml+Lx_down, -Ly-Lpml, 0, lc, 15);
+  gmsh::model::geo::addPoint(Lx+Lpml+Lx_down, Ly+Lpml, 0, lc, 16);
+  gmsh::model::geo::addPoint(-Lx-Lpml, Ly+Lpml, 0, lc, 17);
+
+  gmsh::model::geo::addPoint(-Lx-Lpml, Cy, 0, lc, 49);
+  gmsh::model::geo::addPoint(-Lx, Cy, 0, lc, 50);
+  gmsh::model::geo::addPoint(Lx+Lx_down, Cy, 0, lc, 51);
+  gmsh::model::geo::addPoint(Lx+Lpml+Lx_down, Cy, 0, lc, 52);
+  
+  // adding the lines
+  int l50 = gmsh::model::geo::addLine(49,50);
+  int l51 = gmsh::model::geo::addLine(50,51);
+  int l52 = gmsh::model::geo::addLine(51,52);
+  
+  int l1 = gmsh::model::geo::addLine(2,3);
+  int l2 = gmsh::model::geo::addLine(3,524);
+  int l2222 = gmsh::model::geo::addLine(524,51);
+  int l222 = gmsh::model::geo::addLine(51,4);
+  int l3 = gmsh::model::geo::addLine(4,5);
+  
+  int l4 = gmsh::model::geo::addLine(50,523);
+  int l4444 = gmsh::model::geo::addLine(523,2);
+  int l124 = gmsh::model::geo::addLine(5,50);
+  int l5 = gmsh::model::geo::addLine(7,9);
+  int l6 = gmsh::model::geo::addLine(9,15);
+  int l7 = gmsh::model::geo::addLine(15,8);
+  int l8 = gmsh::model::geo::addLine(8,52);
+  
+  int l888 = gmsh::model::geo::addLine(52,10);
+  int l9 = gmsh::model::geo::addLine(10,16);
+  int l10 = gmsh::model::geo::addLine(16,11);
+  int l11 = gmsh::model::geo::addLine(11,13);
+  int l12 = gmsh::model::geo::addLine(13,17);
+  int l13 = gmsh::model::geo::addLine(17,12);
+
+  int l14 = gmsh::model::geo::addLine(6,49);
+  int l114 = gmsh::model::geo::addLine(49,12);
+  int l15 = gmsh::model::geo::addLine(6,14);
+  int l16 = gmsh::model::geo::addLine(14,7);
+  int l17 = gmsh::model::geo::addLine(3,9);
+  int l18 = gmsh::model::geo::addLine(3,8);
+  int l19 = gmsh::model::geo::addLine(4,10);
+  int l20 = gmsh::model::geo::addLine(4,11);
+
+  int l21 = gmsh::model::geo::addLine(5,13);
+  int l22 = gmsh::model::geo::addLine(5,12);
+  int l23 = gmsh::model::geo::addLine(2,6);
+  int l24 = gmsh::model::geo::addLine(2,7);
+
+  // curve loops
+  int c1 = gmsh::model::geo::addCurveLoop({l1, l2, l2222, -l51, l4, l4444});
+  int p1 = gmsh::model::geo::addPlaneSurface({c1});
+  int c2 = gmsh::model::geo::addCurveLoop({l51, l222, l3, l124});
+  int p2 = gmsh::model::geo::addPlaneSurface({c2});
+
+  // curve loops
+  int c22 = gmsh::model::geo::addCurveLoop({l18, l8, -l52, -l2222, -l2});
+  int p7 = gmsh::model::geo::addPlaneSurface({c22});
+  int c3 = gmsh::model::geo::addCurveLoop({l6, l7, -l18, l17});
+  int p3 = gmsh::model::geo::addPlaneSurface({c3});
+
+  int c4 = gmsh::model::geo::addCurveLoop({l19, l9, l10, -l20});
+  int p4 = gmsh::model::geo::addPlaneSurface({c4});
+  int c5 = gmsh::model::geo::addCurveLoop({l20, l11, -l21, -l3});
+  int p5 = gmsh::model::geo::addPlaneSurface({c5});
+  
+  int c6 = gmsh::model::geo::addCurveLoop({l21, l12, l13, -l22});
+  int p6 = gmsh::model::geo::addPlaneSurface({c6});
+  int c8 = gmsh::model::geo::addCurveLoop({-l24, l16, l15, l23});
+  int p8 = gmsh::model::geo::addPlaneSurface({c8});
+  int c9 = gmsh::model::geo::addCurveLoop({l5, -l17, -l1, l24});
+  int p9 = gmsh::model::geo::addPlaneSurface({c9});
+  
+  int c23 = gmsh::model::geo::addCurveLoop({l52, l888, -l19, -l222});
+  int p10 = gmsh::model::geo::addPlaneSurface({c23});
+  int c24 = gmsh::model::geo::addCurveLoop({-l23, -l14, -l50, -l4, -l4444});
+  int p11 = gmsh::model::geo::addPlaneSurface({c24});
+  int c25 = gmsh::model::geo::addCurveLoop({l50, -l124, l22, -l114});
+  int p12 = gmsh::model::geo::addPlaneSurface({c25});
+
+  // Transfinite curves
+  gmsh::model::geo::synchronize();
+  std::vector<int> lines1{l1, l5, l11, l3};
+  for (unsigned int j=0; j<=lines1.size(); j++) {
+    gmsh::model::geo::mesh::setTransfiniteCurve(lines1[j], Nx, "Progression", 1);
+  }
+  
+  std::vector<int> lines2{l114, l124, l222, l888};
+  for (unsigned int j=0; j<=lines2.size(); j++) {
+    gmsh::model::geo::mesh::setTransfiniteCurve(lines2[j], Ny_short+1, "Progression", 1);
+  }
+  
+  std::vector<int> lines3{l17, l18, l7, l6, l19, l9, l10, l20, l21, l12, l13, l22, l23, l15, l16, l24, l52, l50};
+  for (unsigned int j=0; j<=lines3.size(); j++) {
+    gmsh::model::geo::mesh::setTransfiniteCurve(lines3[j], NrOfLayers+1, "Progression", 1);
+  }
+  
+  
+
+  // embed points
+  gmsh::model::geo::synchronize();
+  gmsh::model::mesh::embed(0, {1,523,524}, 2, p1);
+  for (int j=501; j<=521; j++) {
+    gmsh::model::mesh::embed(0, {j}, 2, p1);
+  }
+
+  // translation
+  std::vector<std::pair<int, int> > out{{2,p1}, {2,p2}, {2,p3}, {2,p4}, {2,p5}, {2,p6}, {2,p7}, {2,p8}, {2,p9}, {2,p10}, {2,p11}, {2,p12}};
+  gmsh::model::geo::translate(out, 2.*sigma, 3.*sigma, 0.);
+  
+  // physical groups
+  gmsh::model::addPhysicalGroup(1, {l1, l2, l2222, l222, l3, l124, l4, l4444}, 1);
+  gmsh::model::setPhysicalName(1, 1, "gamma");
+  gmsh::model::addPhysicalGroup(1, {l51}, 2);
+  gmsh::model::setPhysicalName(1, 2, "control_line");
+
+  gmsh::model::addPhysicalGroup(2, {p1,p2}, 1);
+  gmsh::model::setPhysicalName(2, 1, "omega");
+  gmsh::model::addPhysicalGroup(2, {p3, p4, p5, p6, p7, p8, p9, p10, p11, p12}, 2);
+  gmsh::model::setPhysicalName(2, 2, "omega_pml");
+
+  // generate mesh
+  gmsh::model::mesh::generate(2);
+  gmsh::model::mesh::setOrder(MeshElemOrder);
+  gmsh::write("hotjet.msh");
+
+  // pass physical domain boundary coordinates for gmshfem pml
+  ylim.first = -Cy;
+  ylim.second = Ly+Cy;
+  xlim.first = -5*sigma;
+  xlim.second = 18*sigma;
+  // gmsh::finalize();
+}
\ No newline at end of file
diff --git a/examples/helmholtzFlow/hotjet/mesh.h b/examples/helmholtzFlow/hotjet/mesh.h
new file mode 100644
index 0000000000000000000000000000000000000000..5eeab1ce3157e721d0eb6d3e3f8ef293f76db9e3
--- /dev/null
+++ b/examples/helmholtzFlow/hotjet/mesh.h
@@ -0,0 +1,6 @@
+#include <gmsh.h>
+#include <math.h>
+#include <iostream>
+
+void meshJet(const double sigma, const double SourceMeshRefinement, const int NrOfLayers, const int MeshElemOrder, 
+    std::pair<double, double> &xlim, std::pair<double, double> &ylim, double &lc);
\ No newline at end of file