diff --git a/blast/coupler.py b/blast/coupler.py
index 58430e8e2d4bbaabc7fb114a7cf9d0e7551b0f14..8d826ca0cd8304bba459ae303415647f6cacda8d 100644
--- a/blast/coupler.py
+++ b/blast/coupler.py
@@ -22,6 +22,7 @@ import fwk
 from fwk.coloring import ccolors
 import math
 import numpy as np
+from blast.utils import getSolution
 
 class Coupler:
     def __init__(self, iSolverAPI, vSolver, _maxCouplIter=150, _couplTol=1e-4, _iterPrint=1, _resetInv=False, sfx=''):
@@ -68,6 +69,8 @@ class Coupler:
         cdPrev = 0.0
 
         while couplIter < self.maxIter:
+            #if(couplIter%5 == 0 or couplIter == 0):
+                #print("Inviscid solving !!")
             # Impose blowing boundary condition in the inviscid solver.
             self.tms['processing'].start()
             self.isol.getBlowingBoundary(couplIter)
@@ -76,6 +79,7 @@ class Coupler:
             self.tms['processing'].stop()
 
             # Run inviscid solver.
+            
             self.tms['inviscid'].start()
             if self.resetInviscid:
                 self.isol.reset()
@@ -89,7 +93,8 @@ class Coupler:
                 if write:
                     self.isol.writeCp(sfx='_inviscid'+self.filesfx)
 
-            # Impose inviscid boundary in the viscous solver.
+            #if(couplIter%5 == 0 or couplIter == 0):
+                # Impose inviscid boundary in the viscous solver.
             self.tms['processing'].start()
             self.isol.getInviscidBC()
             self.isol.interpolate('i2v')
@@ -98,8 +103,28 @@ class Coupler:
 
             # Run viscous solver.
             self.tms['viscous'].start()
-            vEc = self.vsol.run()
+            if couplIter == 0:
+                inverse = 1
+            else:
+                inverse = 1
+            
+            vEc = self.vsol.run(inverse,couplIter)
+            inverse_sol = getSolution(self.isol.sec,True,toW ='all',sfx = str(inverse))
             self.tms['viscous'].stop()
+            
+            from matplotlib import pyplot as plt
+            plt.figure()
+            plt.plot(inverse_sol[0]['x'][458:], inverse_sol[0]['deltaStar'][458:])
+            plt.title('deltaStar in the wake')
+            plt.draw()
+            plt.figure()
+            plt.plot(inverse_sol[0]['x'][458:], inverse_sol[0]['ue'][458:],label = "ue")
+            plt.title('Edge velocity in the wake')
+            plt.plot(inverse_sol[0]['x'][458:], inverse_sol[0]['vt'][458:],label = "vt")
+            plt.legend()
+            plt.draw()
+
+            plt.show()
 
             aeroCoeffs['Cl'].append(self.isol.getCl())
             aeroCoeffs['Cd'].append(self.isol.getCd() + self.vsol.Cdf)
diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index b7542bc56ab9b64d9600921f416f5607e88d3053..88b43e4f086d55ea6942b3ca03ce2a8ab452da23 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/blast/interfaces/DSolversInterface.py
@@ -66,6 +66,112 @@ class SolversInterface:
                 reg.setxxExt(self.xxExt[iSec][i])
                 reg.setVtExt(self.vtExt[iSec][i])
 
+        """
+        with open('viscous_input_mach.txt','w') as file:
+            for iSec in range(self.cfg['nSections']):
+                for i, reg in enumerate(self.sec[iSec]):
+                    if reg.name == 'wake':
+                        iReg = 1
+                    elif reg.name == 'lower' or reg.name == 'upper':
+                        iReg = 0
+                    else:
+                        raise RuntimeError('Invalid region', reg.name)
+                    
+                    if reg.name == 'upper':
+                        for MachE in self.vBnd[iSec][iReg].getMach(reg.name)[:]:                  
+                            file.write(f"{MachE}\n")
+
+        with open('boundary_coordinates.txt','w') as file:
+            for iSec in range(self.cfg['nSections']):
+                for i, reg in enumerate(self.sec[iSec]):
+                    if reg.name == 'wake':
+                        iReg = 1
+                    elif reg.name == 'lower' or reg.name == 'upper':
+                        iReg = 0
+                    else:
+                        raise RuntimeError('Invalid region', reg.name)
+                    
+                    if reg.name == 'upper':
+                        nodesCoord = self.vBnd[iSec][0].nodesCoord
+                        for point in nodesCoord:
+                                x, y,_,index = point  # Assuming point is a [x, y] list or tuple
+                                if y >= 0:
+                                    file.write(f"{x} {y} {index}\n")  # Writing x and y in two columns
+
+
+        with open('viscous_input_rho.txt','w') as file:
+            for iSec in range(self.cfg['nSections']):
+                for i, reg in enumerate(self.sec[iSec]):
+                    if reg.name == 'wake':
+                        iReg = 1
+                    elif reg.name == 'lower' or reg.name == 'upper':
+                        iReg = 0
+                    else:
+                        raise RuntimeError('Invalid region', reg.name)
+                    
+                    if reg.name == 'upper':
+                        for RhoE in self.vBnd[iSec][iReg].getRho(reg.name)[:]:
+                  
+                            file.write(f"{RhoE}\n")
+
+        with open('viscous_input_vt.txt','w') as file:
+            for iSec in range(self.cfg['nSections']):
+                for i, reg in enumerate(self.sec[iSec]):
+                    if reg.name == 'wake':
+                        iReg = 1
+                    elif reg.name == 'lower' or reg.name == 'upper':
+                        iReg = 0
+                    else:
+                        raise RuntimeError('Invalid region', reg.name)
+                    
+                    if reg.name == 'upper':
+                        for VtE in self.vBnd[iSec][iReg].getVt(reg.name)[:]:
+                            file.write(f"{VtE}\n")
+
+        with open('viscous_input_deltaStarExt.txt','w') as file:
+            for iSec in range(self.cfg['nSections']):
+                for i, reg in enumerate(self.sec[iSec]):
+                    if reg.name == 'wake':
+                        iReg = 1
+                    elif reg.name == 'lower' or reg.name == 'upper':
+                        iReg = 0
+                    else:
+                        raise RuntimeError('Invalid region', reg.name)
+                    
+                    if reg.name == 'upper':
+                        for deltaStarExt in self.deltaStarExt[iSec][i]:
+                            file.write(f"{deltaStarExt}\n")
+
+
+        with open('viscous_input_xxExt.txt','w') as file:
+            for iSec in range(self.cfg['nSections']):
+                for i, reg in enumerate(self.sec[iSec]):
+                    if reg.name == 'wake':
+                        iReg = 1
+                    elif reg.name == 'lower' or reg.name == 'upper':
+                        iReg = 0
+                    else:
+                        raise RuntimeError('Invalid region', reg.name)
+                    
+                    if reg.name == 'upper':
+                        for xxExt in self.xxExt[iSec][i]:
+                            file.write(f"{xxExt}\n")
+
+        with open('viscous_input_vtExt.txt','w') as file:
+            for iSec in range(self.cfg['nSections']):
+                for i, reg in enumerate(self.sec[iSec]):
+                    if reg.name == 'wake':
+                        iReg = 1
+                    elif reg.name == 'lower' or reg.name == 'upper':
+                        iReg = 0
+                    else:
+                        raise RuntimeError('Invalid region', reg.name)
+                    
+                    if reg.name == 'upper':
+                        for vtExt in self.vtExt[iSec][i]:
+                            file.write(f"{vtExt}\n")
+        """
+
     def getBlowingBoundary(self, couplIter):
         """Get blowing boundary condition from the viscous solver.
         args:
@@ -84,8 +190,8 @@ class SolversInterface:
 
                 for i in range(len(self.sec[iSec])):
                     self.xxExt[iSec][i] = np.asarray(self.sec[iSec][i].loc)
-                    self.deltaStarExt[iSec][i] = np.asarray(self.sec[iSec][i].deltaStar)
-                    self.vtExt[iSec][i] = np.asarray(self.sec[iSec][i].vt)
+                    self.deltaStarExt[iSec][i] = np.asarray(self.sec[iSec][i].deltaStar) #deltaStarExt =
+                    self.vtExt[iSec][i] = np.asarray(self.sec[iSec][i].vt) #.u[3::5] ou vt
     
     def interpolate(self, dir):
         if dir == 'i2v':
@@ -214,7 +320,6 @@ class SolversInterface:
 
                 self.sec[i].append(blast.BoundaryLayer(xtrF, name))
                 self.vSolver.addSection(i, self.sec[i][j])
-
         self.getInviscidBC()
         self.interpolator.inviscidToViscous(self.iBnd, self.vBnd)
 
diff --git a/blast/src/DBoundaryLayer.cpp b/blast/src/DBoundaryLayer.cpp
index 148ed29c5834bfe3e8ba4b65b5e17afddb6e066c..63e0739dc25668c00f49a27736c0196e4ae28ba2 100644
--- a/blast/src/DBoundaryLayer.cpp
+++ b/blast/src/DBoundaryLayer.cpp
@@ -115,23 +115,42 @@ void BoundaryLayer::reset() {
  *
  * @param Re Freestream Reynolds number.
  */
-void BoundaryLayer::setStagnationBC(double Re) {
-  double dv0 = std::abs((vt[1] - vt[0]) / (loc[1] - loc[0]));
-  if (dv0 > 0)
-    u[0] = sqrt(0.075 / (Re * dv0)); // Theta
-  else
-    u[0] = 1e-6; // Theta
-  u[1] = 2.23;   // H
-  u[2] = 0.;     // N
-  u[3] = vt[0];  // ue
-  u[4] = 0.;     // Ctau
-
-  Ret[0] = u[3] * u[0] * Re;
-  if (Ret[0] < 1.) {
-    Ret[0] = 1.;
-    u[3] = Ret[0] / (u[0] * Re);
+void BoundaryLayer::setStagnationBC(double Re, int inverse, int couplIter) {
+    if(inverse == 0 || couplIter == 0){
+    double dv0 = std::abs((vt[1] - vt[0]) / (loc[1] - loc[0]));
+    if (dv0 > 0)
+      u[0] = sqrt(0.075 / (Re * dv0)); // Theta
+    else
+      u[0] = 1e-6; // Theta
+    u[1] = 2.23;   // H
+    u[2] = 0.;     // N
+    u[3] = vt[0];  // ue
+    u[4] = 0.;     // Ctau
+
+    Ret[0] = u[3] * u[0] * Re;
+    if (Ret[0] < 1.) {
+      Ret[0] = 1.;
+      u[3] = Ret[0] / (u[0] * Re);
+    }
+    deltaStar[0] = u[0] * u[1];
+  }
+
+  
+  if(inverse == 1 && couplIter > 0){
+    u[1] = 2.23;   // H
+    u[2] = 0.;     // N
+    u[3] = vt[0];  // ue
+    u[4] = 0.;     // Ctau
+
+    u[0] = deltaStar[0] / u[1]; //Theta
+
+    Ret[0] = u[3] * u[0] * Re;
+    if (Ret[0] < 1.) {
+      Ret[0] = 1.;
+      u[3] = Ret[0] / (u[0] * Re);
+    }  
+  
   }
-  deltaStar[0] = u[0] * u[1];
 
   std::vector<double> lamParam(8, 0.);
   closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], Me[0], name);
@@ -154,25 +173,41 @@ void BoundaryLayer::setStagnationBC(double Re) {
  * @param LowerCond Variables at the last point on the lower side.
  */
 void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond,
-                              std::vector<double> LowerCond) {
+                              std::vector<double> LowerCond, double deltaStarTeUpper, double deltaStarTeLower,int inverse, int couplIter) {
   if (name != "wake")
     std::cout << "Warning: Imposing wake boundary condition on " << name
               << std::endl;
-  u[0] = UpperCond[0] + LowerCond[0]; // (Theta_up+Theta_low).
-  u[1] = (UpperCond[0] * UpperCond[1] + LowerCond[0] * LowerCond[1]) /
-         u[0];  // ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low).
-  u[2] = 9.;    // Amplification factor.
-  u[3] = vt[0]; // Inviscid velocity.
-  u[4] = (UpperCond[0] * UpperCond[4] + LowerCond[0] * LowerCond[4]) /
-         u[0]; // ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low).
 
-  Ret[0] = u[3] * u[0] * Re; // Reynolds number based on the momentum thickness.
+  if(inverse == 0 || couplIter == 0){
+    u[0] = UpperCond[0] + LowerCond[0]; // (Theta_up+Theta_low).
+    u[1] = (UpperCond[0] * UpperCond[1] + LowerCond[0] * LowerCond[1]) /
+          u[0];  // ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low).
+    u[2] = 9.;    // Amplification factor.
+    u[3] = vt[0]; // Inviscid velocity.
+    u[4] = (UpperCond[0] * UpperCond[4] + LowerCond[0] * LowerCond[4]) /
+          u[0]; // ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low).
+
+    deltaStar[0] = u[0] * u[1];
+    Ret[0] = u[3] * u[0] * Re; // Reynolds number based on the momentum thickness.
+  }
+
+
+  if(inverse == 1 && couplIter > 0){
+    deltaStar[0] = deltaStarTeUpper + deltaStarTeLower; // deltaStar_up+deltaStar_low
+    u[0] = UpperCond[0] + LowerCond[0]; // Theta_up+Theta_low
+    u[1] = deltaStar[0] / u[0]; // H = deltaStar/Theta 
+    u[2] = 9.;    // Amplification factor.
+    u[3] = vt[0]; // Inviscid velocity.
+    u[4] = (UpperCond[0] * UpperCond[4] + LowerCond[0] * LowerCond[4]) /
+          u[0]; // ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low).
+    Ret[0] = u[3] * u[0] * Re; // Reynolds number based on the momentum thickness.
+  }
 
   if (Ret[0] < 1.) {
     Ret[0] = 1.;
     u[3] = Ret[0] / (u[0] * Re);
   }
-  deltaStar[0] = u[0] * u[1];
+
 
   // Laminar closures.
   std::vector<double> lamParam(8, 0.);
@@ -190,10 +225,45 @@ void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond,
     regime[k] = 1;
 }
 
+/**
+ * @brief Compute the new etimate of deltaStar in each station of the region, when using semi-inverse method.
+ *
+ */
+
+void BoundaryLayer::updateDeltaStar() {
+
+  for (size_t i = 1; i < nNodes; ++i) {
+    double omega = 0.1; //Relaxation factor
+    double deltaStar_prev = deltaStar[i]; 
+       
+    deltaStar[i] = deltaStar_prev + deltaStar_prev * omega * ((u[i * nVar + 3] / vt[i]) - 1); //Update formula from Carter
+
+
+    if (std::isnan(deltaStar[i - 1])) {
+      if (rhoe[i] == 0.)
+        std::cout << "Density is zero at point " << i << std::endl;
+      if ((loc[i] - loc[i - 1]) == 0.)
+        std::cout << "Points " << i - 1 << " and " << i
+                  << " are at the same position " << loc[i]
+                  << ", resulting in an infinite gradient." << std::endl;
+      if (std::isnan(deltaStar[i]))
+        std::cout << "Displacement thickness is nan at position " << i
+                  << std::endl;
+      if (std::isnan(deltaStar[i - 1]))
+        std::cout << "Displacement thickness is nan at position " << i - 1
+                  << std::endl;
+      throw std::runtime_error(
+          "NaN detected in deltaStar update(see reason above) ");
+
+    }
+  }
+}
+
 /**
  * @brief Compute the blowing velocity in each station of the region.
  *
  */
+
 void BoundaryLayer::computeBlowingVelocity() {
   blowingVelocity.resize(nNodes - 1, 0.);
   for (size_t i = 1; i < nNodes; ++i) {
diff --git a/blast/src/DBoundaryLayer.h b/blast/src/DBoundaryLayer.h
index 87102059ad31837068fbb5396c675b511ff1115a..172117c515b69f514abaafb06c011782f00492d4 100644
--- a/blast/src/DBoundaryLayer.h
+++ b/blast/src/DBoundaryLayer.h
@@ -4,6 +4,7 @@
 #include "DClosures.h"
 #include "blast.h"
 #include "wObject.h"
+#include <vector>
 
 #include <algorithm>
 
@@ -137,9 +138,9 @@ public:
   void setxtrF(double const _xtrF) { xtrF = _xtrF; };
 
   // Boundary conditions.
-  void setStagnationBC(double Re);
+  void setStagnationBC(double Re, int inverse, int couplIter);
   void setWakeBC(double Re, std::vector<double> upperConditions,
-                 std::vector<double> lowerConditions);
+                 std::vector<double> lowerConditions, double deltaStarTeUpper, double deltaStarTeLower,int inverse, int couplIter);
 
   // Getters.
   size_t getnVar() const { return nVar; };
@@ -148,6 +149,7 @@ public:
 
   // Others
   void printSolution(size_t iPoint) const;
+  void updateDeltaStar();
   void computeBlowingVelocity();
   std::vector<std::vector<double>> evalGradwrtRho();
   std::vector<std::vector<double>> evalGradwrtVt();
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index 575682b2ab8b837a4a96295993f853ffdf85e7c4..23b6657c638a9b580d4492e059251df0e36073c4 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -1,6 +1,7 @@
 #include "DDriver.h"
 #include "DBoundaryLayer.h"
 #include "DSolver.h"
+#include "DSolverInverse.h"
 #include <iomanip>
 
 #define ANSI_COLOR_RED "\x1b[1;31m"
@@ -47,6 +48,8 @@ Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
 
   // Time solver initialization.
   tSolver = new Solver(_CFL0, _Minf, Re);
+  tSolverInverse = new SolverInverse(_CFL0, _Minf, Re);
+  std::cout << "SolverInverse Re : " << tSolverInverse->Re << std::endl; 
 
   // Numerical parameters.
   std::cout << "--- Viscous solver setup ---\n"
@@ -62,6 +65,7 @@ Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
  */
 Driver::~Driver() {
   delete tSolver;
+  delete tSolverInverse;
   std::cout << "~blast::Driver()\n";
 } // end ~Driver
 
@@ -79,7 +83,7 @@ void Driver::reset() {
  *
  * @return int Solver exit code.
  */
-int Driver::run() {
+int Driver::run(int inverse, int couplIter) {
   bool lockTrans = false; // Flag activation of transition routines.
   int pointExitCode; // Output of pseudo time integration (0: converged; -1:
                      // unsuccessful, failed; >0: unsuccessful nb iter).
@@ -113,7 +117,7 @@ int Driver::run() {
       if (reg->name == "upper" || reg->name == "lower") {
         for (size_t k = 0; k < reg->nNodes; k++)
           reg->regime[k] = 0;
-        reg->setStagnationBC(Re);
+        reg->setStagnationBC(Re,inverse,couplIter);
       }
 
       else if (reg->name == "wake") // Wake
@@ -134,19 +138,30 @@ int Driver::run() {
                                        sections[iSec][1]->getnVar() +
                                    k];
         }
-        reg->setWakeBC(Re, upperConditions, lowerConditions);
+        double deltaStarTeUpper = sections[iSec][0]->deltaStar[sections[iSec][0]->nNodes - 1];
+        double deltaStarTeLower = sections[iSec][1]->deltaStar[sections[iSec][1]->nNodes - 1]; 
+        reg->setWakeBC(Re, upperConditions, lowerConditions, deltaStarTeUpper, deltaStarTeLower,inverse,couplIter);
         lockTrans = true;
       } else
         throw std::runtime_error("Wrong region name\n");
 
       // Loop over points
       for (size_t iPoint = 1; iPoint < reg->nNodes; ++iPoint) {
-        // Initialize solution at point
-        tSolver->initialCondition(iPoint, reg);
         // Solve equations
-        tms["1-Solver"].start();
-        pointExitCode = tSolver->integration(iPoint, reg);
-        tms["1-Solver"].stop();
+        if(inverse == 0) {
+          tms["1-Solver"].start();
+          tSolver->initialCondition(iPoint, reg); // Initialize solution at point
+          pointExitCode = tSolver->integration(iPoint, reg);
+          tms["1-Solver"].stop();
+        }
+        
+        if(inverse == 1) {
+          tms["1-Solver Inverse"].start();
+          tSolverInverse->initialCondition(iPoint, reg, couplIter); // Initialize solution at point
+          //std::cout << "Inverse mode activated at node : " << iPoint << " of " << reg->name << " region. " << std::endl;
+          pointExitCode = tSolverInverse->solveInverse(iPoint, reg);
+          tms["1-Solver Inverse"].stop();
+        }
 
         if (pointExitCode != 0)
           convergenceStatus[iSec][iRegion].push_back(iPoint);
@@ -157,7 +172,7 @@ int Driver::run() {
           // Forced transition
           if (reg->xtrF != -1 && reg->xoc[iPoint] >= reg->xtrF) {
             reg->u[iPoint * reg->getnVar() + 2] = reg->getnCrit();
-            averageTransition(iPoint, reg, true);
+            averageTransition(iPoint, reg, true, inverse);
             if (verbose > 1) {
               std::cout << std::fixed;
               std::cout << std::setprecision(2);
@@ -172,8 +187,9 @@ int Driver::run() {
 
             // Amplification factor is compared to critical amplification
             // factor.
+
             if (reg->u[iPoint * reg->getnVar() + 2] >= reg->getnCrit()) {
-              averageTransition(iPoint, reg, false);
+              averageTransition(iPoint, reg, false, inverse);
               if (verbose > 1) {
                 std::cout << std::fixed;
                 std::cout << std::setprecision(2);
@@ -185,7 +201,11 @@ int Driver::run() {
         }
         tms["2-Transition"].stop();
       }
+
       tms["3-Blowing"].start();
+      if(inverse == 1){
+        reg->updateDeltaStar();
+      }
       reg->computeBlowingVelocity();
       tms["3-Blowing"].stop();
       iRegion++;
@@ -235,7 +255,7 @@ void Driver::checkWaves(size_t iPoint, BoundaryLayer *bl) {
  * @param bl BoundaryLayer region.
  * @param forced Flag for forced transition.
  */
-void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced) {
+void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced, int inverse) {
   // Averages solution on transition marker.
   size_t nVar = bl->getnVar();
 
@@ -288,14 +308,25 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced) {
   // Avoid starting with ill conditioned IC for turbulent computation. (Regime
   // was switched above). These initial guess values do not influence the
   // converged solution.
+  
   bl->u[iPoint * nVar + 4] =
-      bl->u[(iPoint - 1) * nVar + 4]; // IC of transition point = turbulent BC
-                                      // (imposed @ previous point).
+        bl->u[(iPoint - 1) * nVar + 4]; // IC of transition point = turbulent BC
+                                        // (imposed @ previous point).
   bl->u[iPoint * nVar + 1] =
-      1.515; // Because we expect a significant drop of the shape factor.
+      1.225; // Because we expect a significant drop of the shape factor.
+  
 
   // Solve point in turbulent configuration.
-  int exitCode = tSolver->integration(iPoint, bl);
+  // Need to be adapted for semi inverse
+  int exitCode = -1;
+
+  if(inverse == 1){
+    exitCode = tSolverInverse->solveInverse(iPoint, bl);
+  }
+
+  else if(inverse == 0){
+    exitCode = tSolver->integration(iPoint, bl);
+  }
 
   if (exitCode != 0)
     std::cout << "Warning: Transition point turbulent computation terminated "
@@ -312,9 +343,11 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced) {
       avLam * lamSol[3] + avTurb * bl->u[(iPoint)*nVar + 3];    // ue.
   bl->u[iPoint * nVar + 4] = avTurb * bl->u[(iPoint)*nVar + 4]; // Ct.
 
-  bl->deltaStar[iPoint - 1] =
-      bl->u[(iPoint - 1) * nVar + 0] * bl->u[(iPoint - 1) * nVar + 1];
-  bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * bl->u[iPoint * nVar + 1];
+  //if(inverse == 0){
+    bl->deltaStar[iPoint - 1] =
+        bl->u[(iPoint - 1) * nVar + 0] * bl->u[(iPoint - 1) * nVar + 1];
+    bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * bl->u[iPoint * nVar + 1];
+  //}
 
   // Recompute closures. (The turbulent BC @ iPoint - 1 does not influence
   // laminar closure relations).
@@ -624,4 +657,4 @@ std::vector<std::vector<double>> Driver::getSolution(size_t iSec) {
   if (verbose > 2)
     std::cout << "Solution structure sent." << std::endl;
   return Sln;
-}
\ No newline at end of file
+}
diff --git a/blast/src/DDriver.h b/blast/src/DDriver.h
index 1ff4ba258f00701696a17b595ce0340638a7c41c..1582c07879a2d0ce6252e82f7635b2353db3801f 100644
--- a/blast/src/DDriver.h
+++ b/blast/src/DDriver.h
@@ -30,6 +30,7 @@ private:
   double span; ///< Wing Span (Used only for drag computation, not used if 2D
                ///< case).
   Solver *tSolver;    ///< Pseudo-time solver.
+  SolverInverse *tSolverInverse; ///< Inverse solver 
   int solverExitCode; ///< Exit code of viscous calculations.
   std::vector<std::vector<std::vector<size_t>>>
       convergenceStatus; ///< Vector containing points that did not converged.
@@ -42,7 +43,7 @@ public:
          double xtrF_top = -1., double xtrF_bot = -1., double _span = 0.,
          unsigned int _verbose = 1);
   ~Driver();
-  int run();
+  int run(int inverse = 0,int couplIter = 0);
   void reset();
 
   // Getters.
@@ -59,7 +60,7 @@ public:
 
 private:
   void checkWaves(size_t iPoint, BoundaryLayer *bl);
-  void averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced);
+  void averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced, int inverse);
   void computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i);
   void computeTotalDrag();
   void computeBlwVel();
diff --git a/blast/src/DFluxes.cpp b/blast/src/DFluxes.cpp
index cb9e69d5942fb63c4b57c18f48b9e7099d741f2f..5ecefe22ce37787283dc5c7d2e066300e8671226 100644
--- a/blast/src/DFluxes.cpp
+++ b/blast/src/DFluxes.cpp
@@ -1,6 +1,7 @@
 #include "DFluxes.h"
 #include "DBoundaryLayer.h"
 #include <Eigen/Dense>
+#include "iomanip"
 
 using namespace blast;
 using namespace Eigen;
@@ -49,7 +50,9 @@ MatrixXd Fluxes::computeJacobian(size_t iPoint, BoundaryLayer *bl,
 
 /**
  * @brief Integral boundary layer equation.
- *
+ * ù
+
+mioyv*
  * @param iPoint Marker id.
  * @param bl Boundary layer region.
  * @param u Solution vector at point iPoint.
@@ -84,6 +87,7 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl,
     bl->delta[iPoint] = lamParam[5];
     bl->ctEq[iPoint] = lamParam[6];
     bl->us[iPoint] = lamParam[7];
+
   }
 
   else if (bl->regime[iPoint] == 1) {
@@ -113,7 +117,8 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl,
   double dct_dx = (u[4] - bl->u[(iPoint - 1) * nVar + 4]) / dx;
   double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx;
 
-  double dueExt_dx = (bl->vt[iPoint] - bl->vt[iPoint - 1]) / dxExt;
+  double dueExt_dx = (bl->vt[iPoint] - bl->vt[iPoint - 1]) / dxExt; //due_dksi_old essayer vt_ext 
+
   double ddeltaStarExt_dx =
       (bl->deltaStarExt[iPoint] - bl->deltaStarExt[iPoint - 1]) / dxExt;
 
@@ -123,6 +128,7 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl,
   // Amplification routine.
   double dN_dx = 0.0;
   double ax = 0.0;
+  
   if (bl->xtrF == -1.0 && bl->regime[iPoint] == 0) {
     // No forced transition and laminar regime.
     ax = amplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]);
@@ -278,6 +284,24 @@ VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl,
          (c * u[0] * u[1] + u[3]));
   }
 
+  //Test of the inverse equations
+  // if(iPoint == 1){
+  //   float res0 = 0.;
+  //   float res1 = 0.;
+  //   float res2 = 0.;
+  //   float res3 = 0.;
+    
+  //   res0 = dt_dx + ((2 + u[1] - (Mea*Mea)) * (bl->deltaStar[iPoint]/(u[1]*u[3])) * due_dx) - (cfa/2);
+  //   res1 = ((bl->deltaStar[iPoint]/u[1]) * dHstar_dx) + ((2*Hstar2a + (Hstara * (1-u[1]))) * (bl->deltaStar[iPoint]/(u[1]*u[3])) * due_dx) - (2*cda) + (Hstara*(cfa/2));
+  //   res2 = dN_dx - ax;
+  //   res3 = (2 * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] * dissipRatio) - 2 * deltaa * (4 / (3 * bl->deltaStar[iPoint]) * (cfa / 2 - (((Hka - 1) / (6.7 * Hka * dissipRatio))*((Hka - 1) / (6.7 * Hka * dissipRatio)))) - 1 / u[3] * due_dx);
+
+  //   std::cout << std::fixed << std::setprecision(15) <<"\n" //<< "Res 0 : " << res0 << "\n"
+  //                     //<< "Res 1 : " << res1 << "\n"
+  //                     //<< "Res 2 : " << res2 << "\n"
+  //                     //<< "Res 3 : " << res3  << " u[4] " << u[4] << "\n";
+  //                     //<< "norme : " << sqrt((res0*res0)+(res1*res1)+(res2*res2)+(res3*res3));
+  // }
   return result;
 }
 
diff --git a/blast/src/DFluxesInverse.cpp b/blast/src/DFluxesInverse.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..69ead8d843fe701ce008f8d55777d7eedb1f5fa4
--- /dev/null
+++ b/blast/src/DFluxesInverse.cpp
@@ -0,0 +1,360 @@
+#include "DFluxesInverse.h"
+#include "DBoundaryLayer.h"
+#include <Eigen/Dense>
+#include <vector>
+#include "iomanip" //for print when debbug, to remove
+
+using namespace blast;
+using namespace Eigen;
+
+
+// /**
+//  * @brief Compute residual vector of the integral boundary layer equations.
+//  *
+//  * @param iPoint Marker id.
+//  * @param bl Boundary layer region.
+//  * @return VectorXd
+//  */
+VectorXd FluxesInverse::computeResiduals(size_t iPoint, BoundaryLayer *bl) {
+  size_t nVar = bl->getnVar(); 
+  std::vector<double> u((nVar-1), 0.);
+  for (size_t k = 0; k < u.size(); ++k){ //Since we remove coupling equation
+    u[k] = bl->u[iPoint * nVar + k + 1];
+  }
+  return blLaws(iPoint, bl, u);
+}
+
+/**
+ * @brief Compute Jacobian of the integral boundary layer equations.
+ *
+ * @param iPoint Marker id.
+ * @param bl Boundary layer region.
+ * @param sysRes Residual of the IBL at point iPoint.
+ * @param eps Pertubation from current solution used to compute the Jacobian.
+ * @return MatrixXd
+ */
+MatrixXd FluxesInverse::computeJacobian(size_t iPoint, BoundaryLayer *bl,
+                                 VectorXd const &sysRes, double eps) {
+  size_t nVar = bl->getnVar();
+  MatrixXd JacMatrix = MatrixXd::Zero(4, 4);
+  std::vector<double> uUp(nVar-1, 0.);
+  for (size_t k = 0; k < uUp.size(); ++k) //Since we remove coupling equation, we have one less unknown => strat loop at k = 1 
+    uUp[k] = bl->u[iPoint * nVar + k + 1];
+  
+  double varSave = 0.;
+  for (size_t iVar = 0; iVar < uUp.size(); ++iVar) { 
+    varSave = uUp[iVar];
+    uUp[iVar] += eps;
+    JacMatrix.col(iVar) = (blLaws(iPoint, bl, uUp) - sysRes) / eps;
+    uUp[iVar] = varSave;
+  }
+  return JacMatrix;
+}
+
+
+/**
+ * @brief Integral boundary layer equation.
+ *
+ * @param iPoint Marker id.
+ * @param bl Boundary layer region.
+ * @param u Solution vector at point iPoint.
+ * @return VectorXd
+ */
+VectorXd FluxesInverse::blLaws(size_t iPoint, BoundaryLayer *bl,
+                        std::vector<double> u) {
+  size_t nVar = bl->getnVar(); //Attention = 5 but here we deal with 4 variables 
+
+  // Dissipation ratio.
+  double dissipRatio = 0.;
+  if (bl->name == "wake")
+    dissipRatio = 0.9;
+  else
+    dissipRatio = 1.;
+
+  double theta = 0.;
+  theta = bl->deltaStar[iPoint] / u[0];
+
+  bl->Ret[iPoint] =
+      std::max(u[2] * theta * Re, 1.0); // Reynolds number based on theta.
+
+  // Boundary layer closure relations.
+  if (bl->regime[iPoint] == 0) {
+    // Laminar closure relations.
+    std::vector<double> lamParam(8, 0.);
+    bl->closSolver->laminarClosures(lamParam, theta, u[0], bl->Ret[iPoint],
+                                    bl->Me[iPoint], bl->name);
+    bl->Hstar[iPoint] = lamParam[0];
+    bl->Hstar2[iPoint] = lamParam[1];
+    bl->Hk[iPoint] = lamParam[2];
+    bl->cf[iPoint] = lamParam[3];
+    bl->cd[iPoint] = lamParam[4];
+    bl->delta[iPoint] = lamParam[5];
+    bl->ctEq[iPoint] = lamParam[6];
+    bl->us[iPoint] = lamParam[7];
+  }
+
+  else if (bl->regime[iPoint] == 1) {
+    // Turbulent closure relations.
+    std::vector<double> turbParam(8, 0.);
+    bl->closSolver->turbulentClosures(
+        turbParam, theta, u[0], u[3], bl->Ret[iPoint], bl->Me[iPoint], bl->name);
+    bl->Hstar[iPoint] = turbParam[0];
+    bl->Hstar2[iPoint] = turbParam[1];
+    bl->Hk[iPoint] = turbParam[2];
+    bl->cf[iPoint] = turbParam[3];
+    bl->cd[iPoint] = turbParam[4];
+    bl->delta[iPoint] = turbParam[5];
+    bl->ctEq[iPoint] = turbParam[6];
+    bl->us[iPoint] = turbParam[7];
+  } else {
+    std::cout << "Wrong regime\n" << std::endl;
+    throw;
+  }
+
+  // Gradients computation.
+  double dx = (bl->loc[iPoint] - bl->loc[iPoint - 1]);
+  double due_dx = (u[2] - bl->u[(iPoint - 1) * nVar + 3]) / dx;
+  double dct_dx = (u[3] - (bl->u[(iPoint - 1) * nVar + 4])) / dx;
+  double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx;
+  double dt_dx = ((bl->deltaStar[iPoint]/u[0]) - (bl->deltaStar[iPoint-1]/bl->u[(iPoint - 1)*nVar + 1])) / dx;
+
+
+  // Amplification routine.
+  double dN_dx = 0.0;
+  double ax = 0.0;
+  
+  if (bl->xtrF == -1.0 && bl->regime[iPoint] == 0) {
+    // No forced transition and laminar regime.
+    ax = amplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], theta);
+    dN_dx = (u[1] - bl->u[(iPoint - 1) * nVar + 2]) / dx;
+  }
+
+  double Mea = bl->Me[iPoint];
+  double Hstara = bl->Hstar[iPoint];
+  double Hstar2a = bl->Hstar2[iPoint];
+  double Hka = bl->Hk[iPoint];
+  double cfa = bl->cf[iPoint];
+  double cda = bl->cd[iPoint];
+  double deltaa = bl->delta[iPoint];
+  double ctEqa = bl->ctEq[iPoint];
+
+  //--------------------------------------------------------------------------------//
+  // Integral boundary layer equations. //
+
+  // Space part.
+  // spaceVector(0) = dt_dx + (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx -
+  // cfa / 2; spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 -
+  // u[1])) * u[0] / u[3] * due_dx - 2 * cda + Hstara * cfa / 2; spaceVector(2)
+  // = dN_dx - ax; spaceVector(3) = due_dx - c * (u[1] * dt_dx + u[0] * dH_dx) -
+  // dueExt_dx + cExt * ddeltaStarExt_dx;
+
+  // if (bl->regime[iPoint] == 1)
+  //     spaceVector(4) = (2 * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] *
+  //     dissipRatio) - 2 * deltaa * (4 / (3 * u[1] * u[0]) * (cfa / 2 - (((Hka
+  //     - 1) / (6.7 * Hka * dissipRatio)) * ((Hka - 1) / (6.7 * Hka *
+  //     dissipRatio)))) - 1 / u[3] * due_dx);
+
+  // // Time part.
+  // timeMatrix(0, 0) = u[1] / u[3];
+  // timeMatrix(0, 1) = u[0] / u[3];
+  // timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]);
+
+  // timeMatrix(1, 0) = (1 + u[1] * (1 - Hstara)) / u[3];
+  // timeMatrix(1, 1) = (1 - Hstara) * u[0] / u[3];
+  // timeMatrix(1, 3) = (2 - Hstara * u[1]) * u[0] / (u[3] * u[3]);
+
+  // timeMatrix(2, 2) = 1.;
+
+  // timeMatrix(3, 0) = -c * u[1];
+  // timeMatrix(3, 1) = -c * u[0];
+  // timeMatrix(3, 3) = 1.;
+
+  // if (bl->regime[iPoint] == 1)
+  // {
+  //     timeMatrix(4, 3) = 2 * deltaa / (usa * u[3] * u[3]);
+  //     timeMatrix(4, 4) = 2 * deltaa / (u[3] * u[4] * usa);
+  // }
+  // else
+  //     timeMatrix(4, 4) = 1.0;
+  //--------------------------------------------------------------------------------//
+
+  // VectorXd result = VectorXd::Zero(5);
+
+  // if (bl->regime[iPoint] == 0) {
+  //   result[0] =
+  //       0.5 *
+  //       (-2.0 * u[0] * (1.0 * u[1] - 2.0) *
+  //            (c * (dH_dx * u[0] + dt_dx * u[1]) - cExt * ddeltaStarExt_dx +
+  //             dueExt_dx - due_dx) +
+  //        1.0 * (c * u[0] * u[1] + u[3]) *
+  //            (2.0 * due_dx * u[0] * (2. * Hstar2a - Hstara * (u[1] - 1)) +
+  //             1.0 * u[3] * (Hstara * cfa - 4 * cda + 2. * dHstar_dx * u[0])) +
+  //        1.0 *
+  //            (2. * due_dx * u[0] * (-Mea * Mea + u[1] + 2) +
+  //             u[3] * (-cfa + 2. * dt_dx)) *
+  //            (1.0 * Hstara * c * u[0] * u[1] + 1.0 * Hstara * u[3] -
+  //             2.0 * c * u[0] - 1.0 * u[3])) /
+  //       (c * u[0] * u[1] + u[3]);
+  //   result[1] =
+  //       0.5 *
+  //       (2.0 * u[0] * u[1] * (u[1] - 1) *
+  //            (c * (dH_dx * u[0] + dt_dx * u[1]) - cExt * ddeltaStarExt_dx +
+  //             dueExt_dx - due_dx) -
+  //        1.0 * u[1] * (c * u[0] * u[1] + u[3]) *
+  //            (2. * due_dx * u[0] * (2. * Hstar2a - Hstara * (u[1] - 1)) +
+  //             u[3] * (Hstara * cfa - 4 * cda + 2. * dHstar_dx * u[0])) +
+  //        1.0 *
+  //            (2. * due_dx * u[0] * (-Mea * Mea + u[1] + 2) +
+  //             u[3] * (-cfa + 2. * dt_dx)) *
+  //            (-1.0 * Hstara * c * u[0] * u[1] * u[1] -
+  //             1.0 * Hstara * u[1] * u[3] + 2.0 * c * u[0] * u[1] +
+  //             1.0 * u[1] * u[3] + 1.0 * u[3])) /
+  //       (u[0] * (c * u[0] * u[1] + u[3]));
+  //   result[2] = -1.0 * ax + 1.0 * dN_dx;
+  //   result[3] =
+  //       1.0 * u[3] *
+  //       (-1.0 * c * (dH_dx * u[0] + dt_dx * u[1]) +
+  //        0.5 * c *
+  //            (2. * due_dx * u[0] * (-Mea * Mea + u[1] + 2) +
+  //             u[3] * (-cfa + 2. * dt_dx)) +
+  //        1.0 * cExt * ddeltaStarExt_dx - 1.0 * dueExt_dx + 1.0 * due_dx) /
+  //       (c * u[0] * u[1] + u[3]);
+  //   result[4] = 0.;
+  // } else if (bl->regime[iPoint] == 1) {
+  //   result[0] =
+  //       (-u[0] * (u[1] - 2) *
+  //            (c * (dH_dx * u[0] + dt_dx * u[1]) - cExt * ddeltaStarExt_dx +
+  //             dueExt_dx - due_dx) +
+  //        (c * u[0] * u[1] + u[3]) *
+  //            (2. * due_dx * u[0] * (2. * Hstar2a - Hstara * (u[1] - 1)) +
+  //             u[3] * (Hstara * cfa - 4 * cda + 2. * dHstar_dx * u[0])) /
+  //            2. +
+  //        (2. * due_dx * u[0] * (-Mea * Mea + u[1] + 2) +
+  //         u[3] * (-cfa + 2. * dt_dx)) *
+  //            (Hstara * c * u[0] * u[1] + Hstara * u[3] - 2. * c * u[0] - u[3]) /
+  //            2.) /
+  //       (c * u[0] * u[1] + u[3]);
+  //   result[1] =
+  //       (2. * u[0] * u[1] * (u[1] - 1) *
+  //            (c * (dH_dx * u[0] + dt_dx * u[1]) - cExt * ddeltaStarExt_dx +
+  //             dueExt_dx - due_dx) -
+  //        u[1] * (c * u[0] * u[1] + u[3]) *
+  //            (2. * due_dx * u[0] * (2. * Hstar2a - Hstara * (u[1] - 1)) +
+  //             u[3] * (Hstara * cfa - 4 * cda + 2. * dHstar_dx * u[0])) +
+  //        (2. * due_dx * u[0] * (-Mea * Mea + u[1] + 2) +
+  //         u[3] * (-cfa + 2. * dt_dx)) *
+  //            (-Hstara * c * u[0] * u[1] * u[1] - Hstara * u[1] * u[3] +
+  //             2. * c * u[0] * u[1] + u[1] * u[3] + u[3])) /
+  //       (2. * u[0] * (c 10e-14* u[0] * u[1] + u[3]));
+  //   result[2] = -ax + dN_dx;
+  //   result[3] = u[3] *
+  //               (-2. * c * (dH_dx * u[0] + dt_dx * u[1]) +
+  //                c * (2. * due_dx * u[0] * (-Mea * Mea + u[1] + 2) +
+  //                     u[3] * (-cfa + 2. * dt_dx)) +
+  //                2. * cExt * ddeltaStarExt_dx - 2. * dueExt_dx + 2. * due_dx) /
+  //               (2. * (c * u[0] * u[1] + u[3]));
+  //   result[4] =
+  //       (-Hka * Hka * c * deltaa * dissipRatio * dissipRatio * u[0] * u[1] *
+  //            u[4] *
+  //            (2. * due_dx * u[0] * (-Mea * Mea + u[1] + 2) +
+  //             u[3] * (-cfa + 2. * dt_dx)) /
+  //            2. +
+  //        Hka * Hka * deltaa * dissipRatio * dissipRatio * u[0] * u[1] * u[4] *
+  //            (c * (dH_dx * u[0] + dt_dx * u[1]) - cExt * ddeltaStarExt_dx +
+  //             dueExt_dx - due_dx) +
+  //        usa * (c * u[0] * u[1] + u[3]) *
+  //            (6 * Hka * Hka * dct_dx * deltaa * dissipRatio * dissipRatio *
+  //                 u[0] * u[1] * u[3] +
+  //             16.8 * Hka * Hka * dissipRatio * dissipRatio * u[0] * u[1] *
+  //                 u[3] * u[4] * (-ctEqa + dissipRatio * u[4]) +
+  //             2. * deltaa * u[4] *
+  //                 (3 * Hka * Hka * dissipRatio * dissipRatio * due_dx * u[0] *
+  //                      u[1] +
+  //                  u[3] * (-2. * Hka * Hka * cfa * dissipRatio * dissipRatio +
+  //                          0.0891067052795723 * (Hka - 1) * (Hka - 1)))) /
+  //            6) /
+  //       (Hka * Hka * deltaa * dissipRatio * dissipRatio * u[0] * u[1] *
+  //        (c * u[0] * u[1] + u[3]));
+  // }
+
+  // return result;
+  // -----------------------------------------------------------------------------//
+
+  // --- Inverse system : U = [H,N,ue,c_t] --- //
+
+  VectorXd spaceVector = VectorXd::Zero(4);
+
+
+  spaceVector(0) = dt_dx + ((2 + u[0] - (Mea*Mea)) * (bl->deltaStar[iPoint]/(u[0]*u[2])) * due_dx) - (cfa/2);
+  spaceVector(1) = ((bl->deltaStar[iPoint]/u[0]) * dHstar_dx) + ((2*Hstar2a + (Hstara * (1-u[0]))) * (bl->deltaStar[iPoint]/(u[0]*u[2])) * due_dx) - (2*cda) + (Hstara*(cfa/2));
+  spaceVector(2) = dN_dx - ax;
+
+        
+  if (bl->regime[iPoint] == 1){
+    spaceVector(3) = (2 * deltaa / u[3]) * dct_dx - 5.6 * (ctEqa - u[3] *dissipRatio) - 2 * deltaa * (4 / (3 * u[0] * theta) * (cfa / 2 - (((Hka - 1) / (6.7 * Hka * dissipRatio)) * ((Hka - 1) / (6.7 * Hka * dissipRatio)))) - 1 / u[2] * due_dx);
+  }
+
+  return spaceVector;
+}
+
+/**
+ * @brief Amplification routines of the e^N method for transition capturing.
+ *
+ * @param Hk Kinematic shape parameter.
+ * @param Ret Reynolds number based on the momentum thickness.
+ * @param theta Momentum thickness.
+ * @return double
+ */
+double FluxesInverse::amplificationRoutine(double Hk, double Ret, double theta) const {
+  double dgr = 0.08;
+  double Hk1 = 3.5;
+  double Hk2 = 4.;
+  double Hmi = (1 / (Hk - 1));
+  double logRet_crit =
+      2.492 * std::pow(1. / (Hk - 1.), 0.43) +
+      0.7 * (std::tanh(14. * Hmi - 9.24) +
+             1.0); // value at which the amplification starts to grow
+
+  double ax = 0.;
+  if (Ret <= 0.)
+    Ret = 1.;
+  if (std::log10(Ret) < logRet_crit - dgr)
+    ax = 0.;
+  else {
+    double rNorm = (std::log10(Ret) - (logRet_crit - dgr)) / (2. * dgr);
+    double Rfac = 0.;
+    if (rNorm <= 0)
+      Rfac = 0.;
+    if (rNorm >= 1)
+      Rfac = 1.;
+    else
+      Rfac = 3. * rNorm * rNorm - 2. * rNorm * rNorm * rNorm;
+
+    // Normal envelope amplification rate
+    double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3. * Hmi * Hmi * Hmi +
+               0.1 * std::exp(-20. * Hmi);
+    double arg = 3.87 * Hmi - 2.52;
+    double dndRet = 0.028 * (Hk - 1.) - 0.0345 * std::exp(-(arg * arg));
+    ax = (m * dndRet / theta) * Rfac;
+
+    // Set correction for separated profiles
+    if (Hk > Hk1) {
+      double Hnorm = (Hk - Hk1) / (Hk2 - Hk1);
+      double Hfac = 0.;
+      if (Hnorm >= 1.)
+        Hfac = 1.;
+      else
+        Hfac = 3. * Hnorm * Hnorm - 2. * Hnorm * Hnorm * Hnorm;
+      double ax1 = ax;
+      double gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5.));
+      double tnr = std::tanh(1.2 * (std::log10(Ret) - gro));
+      double ax2 = (0.086 * tnr - 0.25 / std::pow((Hk - 1.), 1.5)) / theta;
+      if (ax2 < 0.)
+        ax2 = 0.;
+      ax = Hfac * ax2 + (1. - Hfac) * ax1;
+      if (ax < 0.)
+        ax = 0.;
+    }
+  }
+  return ax;
+}
diff --git a/blast/src/DFluxesInverse.h b/blast/src/DFluxesInverse.h
new file mode 100644
index 0000000000000000000000000000000000000000..906578d1c973cb13d59cd296ff6fd338ea9acdcd
--- /dev/null
+++ b/blast/src/DFluxesInverse.h
@@ -0,0 +1,32 @@
+#ifndef DFLUXESINVERSE_H
+#define DFLUXESINVERSE_H
+
+#include "DBoundaryLayer.h"
+#include "blast.h"
+#include <Eigen/Dense>
+
+namespace blast {
+
+/**
+ * @brief Residual and Jacobian computation of the unsteady integral boundary
+ * layer equations.
+ * @author Paul Dechamps
+ */
+class BLAST_API FluxesInverse {
+public:
+  double Re; ///< Freestream Reynolds number.
+
+public:
+  FluxesInverse(double _Re) : Re(_Re){};
+  ~FluxesInverse(){};
+  Eigen::VectorXd computeResiduals(size_t iPoint, BoundaryLayer *bl);
+  Eigen::MatrixXd computeJacobian(size_t iPoint, BoundaryLayer *bl,
+                                  Eigen::VectorXd const &sysRes, double eps);
+                         
+private:
+  Eigen::VectorXd blLaws(size_t iPoint, BoundaryLayer *bl,
+                         std::vector<double> u);
+  double amplificationRoutine(double Hk, double Ret, double theta) const;
+};
+} // namespace blast
+#endif // DFLUXESINVERSE_H
diff --git a/blast/src/DSolver.cpp b/blast/src/DSolver.cpp
index 38f46c110ae7acc26979c01e91cf0d2f7a34e644..1228f9ee67ada72dc8c8f46a8c81c06e3c7551f5 100644
--- a/blast/src/DSolver.cpp
+++ b/blast/src/DSolver.cpp
@@ -3,6 +3,7 @@
 #include "DFluxes.h"
 #include <Eigen/Dense>
 #include <iostream>
+#include "iomanip" // To remove , used for debbug
 
 using namespace Eigen;
 using namespace tbox;
@@ -60,6 +61,7 @@ void Solver::initialCondition(size_t iPoint, BoundaryLayer *bl) {
   bl->us[iPoint] = bl->us[iPoint - 1];
   bl->delta[iPoint] = bl->delta[iPoint - 1];
   bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1];
+  
   if (bl->regime[iPoint] == 1 && bl->u[iPoint * nVar + 4] <= 0.)
     bl->u[iPoint * nVar + 4] = 0.03;
 }
@@ -189,9 +191,7 @@ void Solver::resetSolution(size_t iPoint, BoundaryLayer *bl,
   bl->Ret[iPoint] = std::max(bl->u[iPoint * nVar + 3] *
                                  bl->u[iPoint * nVar + 0] * sysMatrix->Re,
                              1.0); // Reynolds number based on theta.
-  bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] *
-                          bl->u[iPoint * nVar + 1]; // Displacement thickness.
-
+  
   if (bl->regime[iPoint] == 0) {
     std::vector<double> lamParam(8, 0.);
     bl->closSolver->laminarClosures(lamParam, bl->u[iPoint * nVar + 0],
diff --git a/blast/src/DSolver.h b/blast/src/DSolver.h
index 604c0061858ee5bcf97404ab4e8f47218ee2781b..aa90b48da73c36d11705dd56c9f822001a1660f3 100644
--- a/blast/src/DSolver.h
+++ b/blast/src/DSolver.h
@@ -24,8 +24,8 @@ private:
   Fluxes *sysMatrix;
 
 public:
-  Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt = 100,
-         double _tol = 1e-8, unsigned int _itFrozenJac = 1);
+  Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt = 500,
+         double _tol = 1e-12, unsigned int _itFrozenJac = 1);
   ~Solver();
   void initialCondition(size_t iPoint, BoundaryLayer *bl);
   int integration(size_t iPoint, BoundaryLayer *bl);
diff --git a/blast/src/DSolverInverse.cpp b/blast/src/DSolverInverse.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4970985949a4b09b3fb0e064b6257cf598b284f8
--- /dev/null
+++ b/blast/src/DSolverInverse.cpp
@@ -0,0 +1,243 @@
+#include "DSolverInverse.h"
+#include "DBoundaryLayer.h"
+#include "DFluxesInverse.h"
+#include <Eigen/Dense>
+#include <iostream>
+#include <iomanip>
+#include <unsupported/Eigen/NonLinearOptimization>
+#include <unsupported/Eigen/NumericalDiff>
+#include <cmath>
+#include <vector>
+
+
+using namespace Eigen;
+using namespace tbox;
+using namespace blast;
+
+/**
+ * @brief Construct a new Time Solver:: Time Solver object.
+ *
+ * @param _CFL0 Initial CFL number.
+ * @param _Minf Freestream Mach number.
+ * @param Re Freestream Reynolds number.
+ * @param _maxIt Maximum number of iterations.
+ * @param _tol Convergence tolerance.
+ * @param _itFrozenJac Number of iterations between which the Jacobian is
+ * frozen.
+ */
+SolverInverse::SolverInverse(double _CFL0, double _Minf, double _Re, unsigned int _maxIt,
+               double _tol, unsigned int _itFrozenJac) {
+  CFL0 = _CFL0;
+  maxIt = _maxIt;
+  tol = _tol;
+  itFrozenJac0 = _itFrozenJac;
+  Minf = std::max(_Minf, 0.1);
+  Re = _Re;
+  sysMatrix = new FluxesInverse(Re);
+}
+
+/**
+ * @brief Destroy the Time Solver:: Time Solver object.
+ *
+ */
+SolverInverse::~SolverInverse() {
+  delete sysMatrix;
+  std::cout << "~blast::SolverInverse()\n";
+}
+
+/**
+ * @brief Impose initial condition at one point.
+ *
+ * @param iPoint Marker id.
+ * @param bl Boundary layer region.
+ */
+void SolverInverse::initialCondition(size_t iPoint, BoundaryLayer *bl, int couplIter) {
+  size_t nVar = bl->getnVar();
+  for (size_t k = 0; k < nVar; ++k)
+    bl->u[iPoint * nVar + k] = bl->u[(iPoint - 1) * nVar + k];
+
+  bl->Ret[iPoint] = bl->Ret[iPoint - 1];
+  bl->cd[iPoint] = bl->cd[iPoint - 1];
+  bl->cf[iPoint] = bl->cf[iPoint - 1];
+  bl->Hstar[iPoint] = bl->Hstar[iPoint - 1];
+  bl->Hstar2[iPoint] = bl->Hstar2[iPoint - 1];
+  bl->Hk[iPoint] = bl->Hk[iPoint - 1];
+  bl->ctEq[iPoint] = bl->ctEq[iPoint - 1];
+  bl->us[iPoint] = bl->us[iPoint - 1];
+  bl->delta[iPoint] = bl->delta[iPoint - 1];
+
+  if(couplIter == 0){
+    if(bl->name == "upper" || bl->name == "lower"){
+      if(bl->xoc[iPoint] == 0){
+        bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1]; 
+      }
+      else{
+        double Rex = Re * bl->xoc[iPoint]; //Self-similarity solution
+        bl->deltaStar[iPoint] = std::max(bl->xoc[iPoint] / sqrt(Rex),bl->deltaStar[iPoint - 1]); //To match the BC in the LE region
+      }
+    }
+    else if(bl->name == "wake"){
+      bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1]; //Currently the only IC that converges for the wake : can be the problem  
+    }
+  }
+  
+  if (bl->regime[iPoint] == 1 && bl->u[iPoint * nVar + 4] <= 0.)
+    bl->u[iPoint * nVar + 4] = 0.03; 
+}
+
+
+/**
+ * @brief Solve steady inverse system with a Newton-Raphson Solver
+ *
+ * @param iPoint Marker id.
+ * @param bl Boundary layer region.
+ * @return int
+ */
+int SolverInverse::solveInverse(size_t iPoint, BoundaryLayer *bl) {
+  size_t nVar = bl->getnVar() ; // /!\ One less unknows since we removed theta from U (no longer coupling equation)
+
+  // Save initial condition.
+  std::vector<double> sln0(nVar, 0.);
+  for (size_t i = 0; i < nVar; ++i)
+    sln0[i] = bl->u[iPoint * nVar + i];
+
+  unsigned int itFrozenJac = itFrozenJac0;
+
+  // Initial system residuals.
+  VectorXd sysRes = sysMatrix->computeResiduals(iPoint, bl);
+
+  double normSysRes0 = sysRes.norm();
+  double normSysRes = normSysRes0;
+  
+  Eigen::MatrixXd JacMatrix = MatrixXd::Zero(4, 4);
+  VectorXd slnIncr = VectorXd::Zero(4);
+
+  unsigned int innerIters = 0; // Inner (non-linear) iterations
+
+  while (innerIters < maxIt) {
+    // Jacobian computation.
+    if (innerIters % itFrozenJac == 0)
+      JacMatrix = sysMatrix->computeJacobian(iPoint, bl, sysRes, 1e-8);
+
+    slnIncr = (JacMatrix).fullPivLu().solve(-sysRes);
+
+    // Solution increment.
+    for (size_t k = 0; k < nVar - 1; ++k){
+      bl->u[iPoint * nVar + k + 1] += slnIncr(k);  
+    }
+
+    bl->u[iPoint * nVar + 1] = std::max(bl->u[iPoint * nVar + 1],  1.0005);
+    bl->u[iPoint * nVar + 0] = std::max((bl->deltaStar[iPoint]/bl->u[iPoint * nVar + 1]), 1e-6);
+
+    sysRes = sysMatrix->computeResiduals(iPoint, bl);
+    normSysRes = sysRes.norm();
+
+    if (normSysRes <= tol * normSysRes0){
+      return 0; // Successfull exit.
+    }
+    
+    innerIters++;
+  }
+
+  if (std::isnan(normSysRes) || normSysRes / normSysRes0 > 1e-3) {
+    resetSolution(iPoint, bl, sln0);
+    return -1;
+  }
+  return innerIters;
+}
+
+
+
+/**
+ * @brief Set time step.
+ *
+ * @param CFL CFL number.
+ * @param dx Cell size.
+ * @param invVel Inviscid velocity.
+ * @return double
+ */
+double SolverInverse::setTimeStep(double CFL, double dx, double invVel) const {
+  // Speed of sound.
+  double gradU2 = (invVel * invVel);
+  double soundSpeed = computeSoundSpeed(gradU2);
+
+  // Velocity of the fastest traveling wave.
+  double advVel = soundSpeed + invVel;
+
+  // Time step computation.
+  return CFL * dx / advVel;
+}
+
+/**
+ * @brief Compute the speed of sound in a cell.
+ *
+ * @param gradU2 Inviscid velocity squared in the cell.
+ * @return double
+ */
+double SolverInverse::computeSoundSpeed(double gradU2) const {
+  return sqrt(1. / (Minf * Minf) + 0.5 * (gamma - 1.) * (1. - gradU2));
+}
+
+/**
+ * @brief Resets the solution of the point wrt its initial condition (sln0) or
+ * the previous point.
+ *
+ * @param iPoint Marker id.
+ * @param bl Boundary layer region.
+ * @param sln0 Initial solution at the point.
+ * @param usePrevPoint if 1, use the solution at previous point instead of sln0.
+ */
+void SolverInverse::resetSolution(size_t iPoint, BoundaryLayer *bl,
+                           std::vector<double> sln0, bool usePrevPoint) {
+  size_t nVar = bl->getnVar();
+
+  // Reset solution.
+  if (usePrevPoint)
+    for (size_t k = 0; k < nVar; ++k)
+      bl->u[iPoint * nVar + k] = bl->u[(iPoint - 1) * nVar + k];
+  else
+    for (size_t k = 0; k < nVar; ++k)
+      bl->u[iPoint * nVar + k] = sln0[k];
+
+  // Reset closures.
+  bl->Ret[iPoint] = std::max(bl->u[iPoint * nVar + 3] *
+                                 bl->u[iPoint * nVar + 0] * sysMatrix->Re,
+                             1.0); // Reynolds number based on theta.
+  //bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] *
+  //                      bl->u[iPoint * nVar + 1]; // Displacement thickness.
+
+  if (bl->regime[iPoint] == 0) {
+    std::vector<double> lamParam(8, 0.);
+    bl->closSolver->laminarClosures(lamParam, bl->u[iPoint * nVar + 0],
+                                    bl->u[iPoint * nVar + 1], bl->Ret[iPoint],
+                                    bl->Me[iPoint], bl->name);
+    bl->Hstar[iPoint] = lamParam[0];
+    bl->Hstar2[iPoint] = lamParam[1];
+    bl->Hk[iPoint] = lamParam[2];
+    bl->cf[iPoint] = lamParam[3];
+    bl->cd[iPoint] = lamParam[4];
+    bl->delta[iPoint] = lamParam[5];
+    bl->ctEq[iPoint] = lamParam[6];
+    bl->us[iPoint] = lamParam[7];
+  } else if (bl->regime[iPoint] == 1) {
+    std::vector<double> turbParam(8, 0.);
+    bl->closSolver->turbulentClosures(
+        turbParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1],
+        bl->u[iPoint * nVar + 4], bl->Ret[iPoint], bl->Me[iPoint], bl->name);
+    bl->Hstar[iPoint] = turbParam[0];
+    bl->Hstar2[iPoint] = turbParam[1];
+    bl->Hk[iPoint] = turbParam[2];
+    bl->cf[iPoint] = turbParam[3];
+    bl->cd[iPoint] = turbParam[4];
+    bl->delta[iPoint] = turbParam[5];
+    bl->ctEq[iPoint] = turbParam[6];
+    bl->us[iPoint] = turbParam[7];
+  }
+}
+
+void SolverInverse::setCFL0(double _CFL0) {
+  if (_CFL0 > 0)
+    CFL0 = _CFL0;
+  else
+    throw std::runtime_error("Negative CFL numbers are not allowed.\n");
+}
diff --git a/blast/src/DSolverInverse.h b/blast/src/DSolverInverse.h
new file mode 100644
index 0000000000000000000000000000000000000000..4b8ce8d636a807c4e76d51403276d123b0ace493
--- /dev/null
+++ b/blast/src/DSolverInverse.h
@@ -0,0 +1,57 @@
+#ifndef DSOLVERINVERSE_H
+#define DSOLVERINVERSE_H
+
+#include "DFluxesInverse.h"
+#include "blast.h"
+#include <Eigen/Dense>
+#include <unsupported/Eigen/NonLinearOptimization>
+
+namespace blast {
+
+/**
+ * @brief Pseudo-time integration.
+ * @author Paul Dechamps
+ */
+class BLAST_API SolverInverse {
+
+private:
+  double CFL0;               ///< Initial CFL number.
+  unsigned int maxIt;        ///< Maximum number of iterations.
+  double tol;                ///< Convergence tolerance.
+  unsigned int itFrozenJac0; ///< Number of iterations between which the
+                             ///< Jacobian is frozen.
+  double const gamma = 1.4;  ///< Air heat capacity ratio.
+  double Minf;               ///< Freestream Mach number.
+
+  FluxesInverse *sysMatrix;
+
+public:
+  SolverInverse(double _CFL0, double _Minf, double _Re, unsigned int _maxIt = 1000,
+         double _tol = 1e-9, unsigned int _itFrozenJac = 1);
+  ~SolverInverse();
+  void initialCondition(size_t iPoint, BoundaryLayer *bl, int couplIter);
+  int solveInverse(size_t iPoint, BoundaryLayer *bl); 
+
+  double Re;               ///< Freestream Mach number.
+
+  // Getters.
+  double getCFL0() const { return CFL0; };
+  unsigned int getmaxIt() const { return maxIt; };
+  unsigned int getitFrozenJac() const { return itFrozenJac0; };
+
+  // Setters.
+  void setCFL0(double _CFL0);
+  void setmaxIt(unsigned int _maxIt) { maxIt = _maxIt; };
+  void setitFrozenJac(unsigned int _itFrozenJac) {
+    itFrozenJac0 = _itFrozenJac;
+  };
+
+private:
+  double setTimeStep(double CFL, double dx, double invVel) const;
+  double computeSoundSpeed(double gradU2) const;
+  void resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double> Sln0,
+                     bool usePrevPoint = false);
+
+};
+} // namespace blast
+#endif // DSOLVERINVERSE_H
diff --git a/blast/src/blast.h b/blast/src/blast.h
index 9b526a72d2ff83cca967adc14a6346f2e42452ae..9aead0e78eebc10bfec9731e8d3458b9f8f1cd61 100644
--- a/blast/src/blast.h
+++ b/blast/src/blast.h
@@ -39,6 +39,7 @@ namespace blast {
 // Solvers
 class Driver;
 class Solver;
+class SolverInverse;
 
 // Adjoint
 class CoupledAdjoint;
diff --git a/blast/thesis_tests/NACA0012.geo b/blast/thesis_tests/NACA0012.geo
new file mode 100644
index 0000000000000000000000000000000000000000..45f9b857e8b81b2865f044f89530e622b48cc111
--- /dev/null
+++ b/blast/thesis_tests/NACA0012.geo
@@ -0,0 +1,389 @@
+/* Airfoil NACA0012 */
+// This file was generated automatically using geoFoils v1.0 (https://github.com/Paul-Dech/geoFoils)
+// @date: 2024-10-09 17:12:32
+// @author: Paul Dechamps
+
+// Geometry
+DefineConstant[ xLgt = { 50.0, Name "Domain length (x-dir)" }  ];
+DefineConstant[ yLgt = { 50.0, Name "Domain length (y-dir)" }  ];
+
+// Mesh
+DefineConstant[ growthRatio = { 1.5, Name "Growth Ratio" }  ];
+DefineConstant[ msTe = { 0.01, Name "Airfoil TE mesh size" }  ];
+DefineConstant[ msLe = { 0.001, Name "Airfoil LE mesh size" }  ];
+
+// Rotation
+DefineConstant[ xRot = { 0.25, Name "Center of rotation" }  ];
+DefineConstant[ angle = { 0*Pi/180, Name "Angle of rotation" }  ];
+Geometry.AutoCoherence = 0; // Needed so that gmsh does not remove duplicate
+
+If (growthRatio == 1.0)
+   msF = msLe;
+Else
+   n = Log(1 - (1 - growthRatio) * xLgt / msLe) / Log(growthRatio);
+   msF = msLe * growthRatio^(n - 1);
+EndIf
+
+/**************
+ Geometry
+ **************/
+Te = 1; // trailing edge
+Le = 151; // leading edge
+
+Point( 1 ) = { 1.0, -0.0, 0.0, msTe };
+Point( 2 ) = { 0.99989034173742, 1.566142976e-05, 0.0};
+Point( 3 ) = { 0.99956141504943, 6.262613794e-05, 0.0};
+Point( 4 ) = { 0.99901336421414, 0.00014083544147, 0.0};
+Point( 5 ) = { 0.99824642962475, 0.0002501917354, 0.0};
+Point( 6 ) = { 0.99726094768414, 0.00039055879171, 0.0};
+Point( 7 ) = { 0.99605735065724, 0.00056176217462, 0.0};
+Point( 8 ) = { 0.99463616648149, 0.00076358976946, 0.0};
+Point( 9 ) = { 0.99299801853525, 0.00099579242133, 0.0};
+Point( 10 ) = { 0.99114362536434, 0.00125808467902, 0.0};
+Point( 11 ) = { 0.9890738003669, 0.0015501456389, 0.0};
+Point( 12 ) = { 0.98678945143658, 0.00187161988278, 0.0};
+Point( 13 ) = { 0.98429158056432, 0.0022221185031, 0.0};
+Point( 14 ) = { 0.98158128339883, 0.00260122020822, 0.0};
+Point( 15 ) = { 0.97865974876603, 0.00300847250009, 0.0};
+Point( 16 ) = { 0.97552825814758, 0.00344339291591, 0.0};
+Point( 17 ) = { 0.97218818511874, 0.0039054703252, 0.0};
+Point( 18 ) = { 0.96864099474595, 0.00439416627328, 0.0};
+Point( 19 ) = { 0.96488824294413, 0.00490891636171, 0.0};
+Point( 20 ) = { 0.96093157579425, 0.00544913165641, 0.0};
+Point( 21 ) = { 0.9567727288213, 0.00601420011357, 0.0};
+Point( 22 ) = { 0.95241352623301, 0.00660348801398, 0.0};
+Point( 23 ) = { 0.94785588011971, 0.00721634139577, 0.0};
+Point( 24 ) = { 0.94310178961561, 0.00785208747637, 0.0};
+Point( 25 ) = { 0.93815334002193, 0.00851003605407, 0.0};
+Point( 26 ) = { 0.93301270189222, 0.00918948088016, 0.0};
+Point( 27 ) = { 0.92768213008025, 0.00988970099293, 0.0};
+Point( 28 ) = { 0.92216396275101, 0.01060996200508, 0.0};
+Point( 29 ) = { 0.91646062035505, 0.01134951733663, 0.0};
+Point( 30 ) = { 0.91057460456685, 0.01210760938604, 0.0};
+Point( 31 ) = { 0.90450849718747, 0.01288347063272, 0.0};
+Point( 32 ) = { 0.8982649590121, 0.01367632466478, 0.0};
+Point( 33 ) = { 0.89184672866292, 0.01448538712671, 0.0};
+Point( 34 ) = { 0.88525662138789, 0.01530986658219, 0.0};
+Point( 35 ) = { 0.87849752782588, 0.01614896528828, 0.0};
+Point( 36 ) = { 0.8715724127387, 0.01700187987791, 0.0};
+Point( 37 ) = { 0.86448431371071, 0.01786780194847, 0.0};
+Point( 38 ) = { 0.8572363398164, 0.01874591855516, 0.0};
+Point( 39 ) = { 0.84983167025668, 0.01963541260873, 0.0};
+Point( 40 ) = { 0.84227355296434, 0.02053546317813, 0.0};
+Point( 41 ) = { 0.83456530317943, 0.02144524569925, 0.0};
+Point( 42 ) = { 0.82671030199505, 0.02236393209232, 0.0};
+Point( 43 ) = { 0.81871199487434, 0.02329069079083, 0.0};
+Point( 44 ) = { 0.81057389013916, 0.02422468668617, 0.0};
+Point( 45 ) = { 0.80229955743119, 0.02516508099279, 0.0};
+Point( 46 ) = { 0.79389262614624, 0.02611103103947, 0.0};
+Point( 47 ) = { 0.78535678384222, 0.02706168999315, 0.0};
+Point( 48 ) = { 0.77669577462167, 0.02801620652246, 0.0};
+Point( 49 ) = { 0.7679133974895, 0.02897372440866, 0.0};
+Point( 50 ) = { 0.75901350468657, 0.02993338211239, 0.0};
+Point( 51 ) = { 0.75, 0.03089431230516, 0.0};
+Point( 52 ) = { 0.74087683705086, 0.03185564137496, 0.0};
+Point( 53 ) = { 0.73164801755993, 0.03281648891572, 0.0};
+Point( 54 ) = { 0.72231758959246, 0.03377596721082, 0.0};
+Point( 55 ) = { 0.71288964578254, 0.03473318072085, 0.0};
+Point( 56 ) = { 0.7033683215379, 0.03568722558628, 0.0};
+Point( 57 ) = { 0.69375779322605, 0.03663718915547, 0.0};
+Point( 58 ) = { 0.68406227634234, 0.03758214954874, 0.0};
+Point( 59 ) = { 0.67428602366091, 0.03852117526877, 0.0};
+Point( 60 ) = { 0.66443332336929, 0.03945332486788, 0.0};
+Point( 61 ) = { 0.65450849718747, 0.0403776466819, 0.0};
+Point( 62 ) = { 0.64451589847224, 0.04129317864062, 0.0};
+Point( 63 ) = { 0.63445991030763, 0.04219894816377, 0.0};
+Point( 64 ) = { 0.62434494358243, 0.04309397215141, 0.0};
+Point( 65 ) = { 0.61417543505533, 0.04397725707679, 0.0};
+Point( 66 ) = { 0.60395584540888, 0.04484779918912, 0.0};
+Point( 67 ) = { 0.59369065729286, 0.04570458483296, 0.0};
+Point( 68 ) = { 0.58338437335805, 0.04654659089014, 0.0};
+Point( 69 ) = { 0.57304151428121, 0.04737278534922, 0.0};
+Point( 70 ) = { 0.56266661678215, 0.04818212800662, 0.0};
+Point( 71 ) = { 0.55226423163383, 0.04897357130255, 0.0};
+Point( 72 ) = { 0.54183892166616, 0.0497460612939, 0.0};
+Point( 73 ) = { 0.53139525976466, 0.05049853876514, 0.0};
+Point( 74 ) = { 0.5209378268646, 0.05122994047724, 0.0};
+Point( 75 ) = { 0.51047120994168, 0.05193920055357, 0.0};
+Point( 76 ) = { 0.5, 0.05262525200057, 0.0};
+Point( 77 ) = { 0.48952879005832, 0.05328702835989, 0.0};
+Point( 78 ) = { 0.4790621731354, 0.05392346548752, 0.0};
+Point( 79 ) = { 0.46860474023534, 0.05453350345464, 0.0};
+Point( 80 ) = { 0.45816107833384, 0.05511608856324, 0.0};
+Point( 81 ) = { 0.44773576836617, 0.05567017546926, 0.0};
+Point( 82 ) = { 0.43733338321785, 0.05619472940438, 0.0};
+Point( 83 ) = { 0.42695848571879, 0.05668872848681, 0.0};
+Point( 84 ) = { 0.41661562664195, 0.05715116611079, 0.0};
+Point( 85 ) = { 0.40630934270714, 0.05758105340305, 0.0};
+Point( 86 ) = { 0.39604415459112, 0.05797742173437, 0.0};
+Point( 87 ) = { 0.38582456494467, 0.05833932527305, 0.0};
+Point( 88 ) = { 0.37565505641757, 0.05866584356678, 0.0};
+Point( 89 ) = { 0.36554008969237, 0.05895608413888, 0.0};
+Point( 90 ) = { 0.35548410152776, 0.05920918508399, 0.0};
+Point( 91 ) = { 0.34549150281253, 0.05942431764848, 0.0};
+Point( 92 ) = { 0.33556667663071, 0.05960068878005, 0.0};
+Point( 93 ) = { 0.32571397633909, 0.05973754363096, 0.0};
+Point( 94 ) = { 0.31593772365766, 0.0598341679995, 0.0};
+Point( 95 ) = { 0.30624220677395, 0.05988989069371, 0.0};
+Point( 96 ) = { 0.2966316784621, 0.05990408580221, 0.0};
+Point( 97 ) = { 0.28711035421746, 0.05987617485653, 0.0};
+Point( 98 ) = { 0.27768241040754, 0.05980562887012, 0.0};
+Point( 99 ) = { 0.26835198244007, 0.05969197023942, 0.0};
+Point( 100 ) = { 0.25912316294914, 0.05953477449294, 0.0};
+Point( 101 ) = { 0.25, 0.059333671875, 0.0};
+Point( 102 ) = { 0.24098649531343, 0.05908834875129, 0.0};
+Point( 103 ) = { 0.2320866025105, 0.05879854882451, 0.0};
+Point( 104 ) = { 0.22330422537833, 0.058464074149, 0.0};
+Point( 105 ) = { 0.21464321615778, 0.05808478593443, 0.0};
+Point( 106 ) = { 0.20610737385376, 0.05766060512952, 0.0};
+Point( 107 ) = { 0.19770044256881, 0.05719151277812, 0.0};
+Point( 108 ) = { 0.18942610986084, 0.05667755014088, 0.0};
+Point( 109 ) = { 0.18128800512566, 0.0561188185773, 0.0};
+Point( 110 ) = { 0.17328969800495, 0.05551547918403, 0.0};
+Point( 111 ) = { 0.16543469682057, 0.05486775218683, 0.0};
+Point( 112 ) = { 0.15772644703566, 0.05417591608471, 0.0};
+Point( 113 ) = { 0.15016832974332, 0.05344030654663, 0.0};
+Point( 114 ) = { 0.1427636601836, 0.05266131506196, 0.0};
+Point( 115 ) = { 0.13551568628929, 0.05183938734803, 0.0};
+Point( 116 ) = { 0.1284275872613, 0.05097502151884, 0.0};
+Point( 117 ) = { 0.12150247217412, 0.05006876602101, 0.0};
+Point( 118 ) = { 0.11474337861211, 0.04912121734409, 0.0};
+Point( 119 ) = { 0.10815327133708, 0.04813301751383, 0.0};
+Point( 120 ) = { 0.1017350409879, 0.04710485137844, 0.0};
+Point( 121 ) = { 0.09549150281253, 0.046037443699, 0.0};
+Point( 122 ) = { 0.08942539543315, 0.04493155605656, 0.0};
+Point( 123 ) = { 0.08353937964495, 0.04378798358951, 0.0};
+Point( 124 ) = { 0.07783603724899, 0.04260755157605, 0.0};
+Point( 125 ) = { 0.07231786991975, 0.04139111187757, 0.0};
+Point( 126 ) = { 0.06698729810778, 0.04013953925954, 0.0};
+Point( 127 ) = { 0.06184665997807, 0.03885372760755, 0.0};
+Point( 128 ) = { 0.05689821038439, 0.03753458605692, 0.0};
+Point( 129 ) = { 0.05214411988029, 0.03618303505469, 0.0};
+Point( 130 ) = { 0.04758647376699, 0.0348000023736, 0.0};
+Point( 131 ) = { 0.0432272711787, 0.03338641909797, 0.0};
+Point( 132 ) = { 0.03906842420575, 0.03194321560181, 0.0};
+Point( 133 ) = { 0.03511175705587, 0.03047131753953, 0.0};
+Point( 134 ) = { 0.03135900525405, 0.0289716418698, 0.0};
+Point( 135 ) = { 0.02781181488126, 0.02744509293315, 0.0};
+Point( 136 ) = { 0.02447174185242, 0.0258925586035, 0.0};
+Point( 137 ) = { 0.02134025123397, 0.02431490653389, 0.0};
+Point( 138 ) = { 0.01841871660117, 0.02271298051574, 0.0};
+Point( 139 ) = { 0.01570841943568, 0.02108759697117, 0.0};
+Point( 140 ) = { 0.01321054856342, 0.01943954159657, 0.0};
+Point( 141 ) = { 0.0109261996331, 0.01776956617535, 0.0};
+Point( 142 ) = { 0.00885637463566, 0.01607838557678, 0.0};
+Point( 143 ) = { 0.00700198146475, 0.01436667495681, 0.0};
+Point( 144 ) = { 0.00536383351851, 0.01263506717596, 0.0};
+Point( 145 ) = { 0.00394264934276, 0.01088415044781, 0.0};
+Point( 146 ) = { 0.00273905231586, 0.009114466231, 0.0};
+Point( 147 ) = { 0.00175357037525, 0.00732650737576, 0.0};
+Point( 148 ) = { 0.00098663578586, 0.00552071653494, 0.0};
+Point( 149 ) = { 0.00043858495057, 0.0036974848482, 0.0};
+Point( 150 ) = { 0.00010965826258, 0.00185715090611, 0.0};
+Point( 151 ) = { 0.0, 0.0, 0.0, msLe };
+Point( 152 ) = { 0.00010965826258, -0.00185715090611, 0.0};
+Point( 153 ) = { 0.00043858495057, -0.0036974848482, 0.0};
+Point( 154 ) = { 0.00098663578586, -0.00552071653494, 0.0};
+Point( 155 ) = { 0.00175357037525, -0.00732650737576, 0.0};
+Point( 156 ) = { 0.00273905231586, -0.009114466231, 0.0};
+Point( 157 ) = { 0.00394264934276, -0.01088415044781, 0.0};
+Point( 158 ) = { 0.00536383351851, -0.01263506717596, 0.0};
+Point( 159 ) = { 0.00700198146475, -0.01436667495681, 0.0};
+Point( 160 ) = { 0.00885637463566, -0.01607838557678, 0.0};
+Point( 161 ) = { 0.0109261996331, -0.01776956617535, 0.0};
+Point( 162 ) = { 0.01321054856342, -0.01943954159657, 0.0};
+Point( 163 ) = { 0.01570841943568, -0.02108759697117, 0.0};
+Point( 164 ) = { 0.01841871660117, -0.02271298051574, 0.0};
+Point( 165 ) = { 0.02134025123397, -0.02431490653389, 0.0};
+Point( 166 ) = { 0.02447174185242, -0.0258925586035, 0.0};
+Point( 167 ) = { 0.02781181488126, -0.02744509293315, 0.0};
+Point( 168 ) = { 0.03135900525405, -0.0289716418698, 0.0};
+Point( 169 ) = { 0.03511175705587, -0.03047131753953, 0.0};
+Point( 170 ) = { 0.03906842420575, -0.03194321560181, 0.0};
+Point( 171 ) = { 0.0432272711787, -0.03338641909797, 0.0};
+Point( 172 ) = { 0.04758647376699, -0.0348000023736, 0.0};
+Point( 173 ) = { 0.05214411988029, -0.03618303505469, 0.0};
+Point( 174 ) = { 0.05689821038439, -0.03753458605692, 0.0};
+Point( 175 ) = { 0.06184665997807, -0.03885372760755, 0.0};
+Point( 176 ) = { 0.06698729810778, -0.04013953925954, 0.0};
+Point( 177 ) = { 0.07231786991975, -0.04139111187757, 0.0};
+Point( 178 ) = { 0.07783603724899, -0.04260755157605, 0.0};
+Point( 179 ) = { 0.08353937964495, -0.04378798358951, 0.0};
+Point( 180 ) = { 0.08942539543315, -0.04493155605656, 0.0};
+Point( 181 ) = { 0.09549150281253, -0.046037443699, 0.0};
+Point( 182 ) = { 0.1017350409879, -0.04710485137844, 0.0};
+Point( 183 ) = { 0.10815327133708, -0.04813301751383, 0.0};
+Point( 184 ) = { 0.11474337861211, -0.04912121734409, 0.0};
+Point( 185 ) = { 0.12150247217412, -0.05006876602101, 0.0};
+Point( 186 ) = { 0.1284275872613, -0.05097502151884, 0.0};
+Point( 187 ) = { 0.13551568628929, -0.05183938734803, 0.0};
+Point( 188 ) = { 0.1427636601836, -0.05266131506196, 0.0};
+Point( 189 ) = { 0.15016832974332, -0.05344030654663, 0.0};
+Point( 190 ) = { 0.15772644703566, -0.05417591608471, 0.0};
+Point( 191 ) = { 0.16543469682057, -0.05486775218683, 0.0};
+Point( 192 ) = { 0.17328969800495, -0.05551547918403, 0.0};
+Point( 193 ) = { 0.18128800512566, -0.0561188185773, 0.0};
+Point( 194 ) = { 0.18942610986084, -0.05667755014088, 0.0};
+Point( 195 ) = { 0.19770044256881, -0.05719151277812, 0.0};
+Point( 196 ) = { 0.20610737385376, -0.05766060512952, 0.0};
+Point( 197 ) = { 0.21464321615778, -0.05808478593443, 0.0};
+Point( 198 ) = { 0.22330422537833, -0.058464074149, 0.0};
+Point( 199 ) = { 0.2320866025105, -0.05879854882451, 0.0};
+Point( 200 ) = { 0.24098649531343, -0.05908834875129, 0.0};
+Point( 201 ) = { 0.25, -0.059333671875, 0.0};
+Point( 202 ) = { 0.25912316294914, -0.05953477449294, 0.0};
+Point( 203 ) = { 0.26835198244007, -0.05969197023942, 0.0};
+Point( 204 ) = { 0.27768241040754, -0.05980562887012, 0.0};
+Point( 205 ) = { 0.28711035421746, -0.05987617485653, 0.0};
+Point( 206 ) = { 0.2966316784621, -0.05990408580221, 0.0};
+Point( 207 ) = { 0.30624220677395, -0.05988989069371, 0.0};
+Point( 208 ) = { 0.31593772365766, -0.0598341679995, 0.0};
+Point( 209 ) = { 0.32571397633909, -0.05973754363096, 0.0};
+Point( 210 ) = { 0.33556667663071, -0.05960068878005, 0.0};
+Point( 211 ) = { 0.34549150281253, -0.05942431764848, 0.0};
+Point( 212 ) = { 0.35548410152776, -0.05920918508399, 0.0};
+Point( 213 ) = { 0.36554008969237, -0.05895608413888, 0.0};
+Point( 214 ) = { 0.37565505641757, -0.05866584356678, 0.0};
+Point( 215 ) = { 0.38582456494467, -0.05833932527305, 0.0};
+Point( 216 ) = { 0.39604415459112, -0.05797742173437, 0.0};
+Point( 217 ) = { 0.40630934270714, -0.05758105340305, 0.0};
+Point( 218 ) = { 0.41661562664195, -0.05715116611079, 0.0};
+Point( 219 ) = { 0.42695848571879, -0.05668872848681, 0.0};
+Point( 220 ) = { 0.43733338321785, -0.05619472940438, 0.0};
+Point( 221 ) = { 0.44773576836617, -0.05567017546926, 0.0};
+Point( 222 ) = { 0.45816107833384, -0.05511608856324, 0.0};
+Point( 223 ) = { 0.46860474023534, -0.05453350345464, 0.0};
+Point( 224 ) = { 0.4790621731354, -0.05392346548752, 0.0};
+Point( 225 ) = { 0.48952879005832, -0.05328702835989, 0.0};
+Point( 226 ) = { 0.5, -0.05262525200057, 0.0};
+Point( 227 ) = { 0.51047120994168, -0.05193920055357, 0.0};
+Point( 228 ) = { 0.5209378268646, -0.05122994047724, 0.0};
+Point( 229 ) = { 0.53139525976466, -0.05049853876514, 0.0};
+Point( 230 ) = { 0.54183892166616, -0.0497460612939, 0.0};
+Point( 231 ) = { 0.55226423163383, -0.04897357130255, 0.0};
+Point( 232 ) = { 0.56266661678215, -0.04818212800662, 0.0};
+Point( 233 ) = { 0.57304151428121, -0.04737278534922, 0.0};
+Point( 234 ) = { 0.58338437335805, -0.04654659089014, 0.0};
+Point( 235 ) = { 0.59369065729286, -0.04570458483296, 0.0};
+Point( 236 ) = { 0.60395584540888, -0.04484779918912, 0.0};
+Point( 237 ) = { 0.61417543505533, -0.04397725707679, 0.0};
+Point( 238 ) = { 0.62434494358243, -0.04309397215141, 0.0};
+Point( 239 ) = { 0.63445991030763, -0.04219894816377, 0.0};
+Point( 240 ) = { 0.64451589847224, -0.04129317864062, 0.0};
+Point( 241 ) = { 0.65450849718747, -0.0403776466819, 0.0};
+Point( 242 ) = { 0.66443332336929, -0.03945332486788, 0.0};
+Point( 243 ) = { 0.67428602366091, -0.03852117526877, 0.0};
+Point( 244 ) = { 0.68406227634234, -0.03758214954874, 0.0};
+Point( 245 ) = { 0.69375779322605, -0.03663718915547, 0.0};
+Point( 246 ) = { 0.7033683215379, -0.03568722558628, 0.0};
+Point( 247 ) = { 0.71288964578254, -0.03473318072085, 0.0};
+Point( 248 ) = { 0.72231758959246, -0.03377596721082, 0.0};
+Point( 249 ) = { 0.73164801755993, -0.03281648891572, 0.0};
+Point( 250 ) = { 0.74087683705086, -0.03185564137496, 0.0};
+Point( 251 ) = { 0.75, -0.03089431230516, 0.0};
+Point( 252 ) = { 0.75901350468657, -0.02993338211239, 0.0};
+Point( 253 ) = { 0.7679133974895, -0.02897372440866, 0.0};
+Point( 254 ) = { 0.77669577462167, -0.02801620652246, 0.0};
+Point( 255 ) = { 0.78535678384222, -0.02706168999315, 0.0};
+Point( 256 ) = { 0.79389262614624, -0.02611103103947, 0.0};
+Point( 257 ) = { 0.80229955743119, -0.02516508099279, 0.0};
+Point( 258 ) = { 0.81057389013916, -0.02422468668617, 0.0};
+Point( 259 ) = { 0.81871199487434, -0.02329069079083, 0.0};
+Point( 260 ) = { 0.82671030199505, -0.02236393209232, 0.0};
+Point( 261 ) = { 0.83456530317943, -0.02144524569925, 0.0};
+Point( 262 ) = { 0.84227355296434, -0.02053546317813, 0.0};
+Point( 263 ) = { 0.84983167025668, -0.01963541260873, 0.0};
+Point( 264 ) = { 0.8572363398164, -0.01874591855516, 0.0};
+Point( 265 ) = { 0.86448431371071, -0.01786780194847, 0.0};
+Point( 266 ) = { 0.8715724127387, -0.01700187987791, 0.0};
+Point( 267 ) = { 0.87849752782588, -0.01614896528828, 0.0};
+Point( 268 ) = { 0.88525662138789, -0.01530986658219, 0.0};
+Point( 269 ) = { 0.89184672866292, -0.01448538712671, 0.0};
+Point( 270 ) = { 0.8982649590121, -0.01367632466478, 0.0};
+Point( 271 ) = { 0.90450849718747, -0.01288347063272, 0.0};
+Point( 272 ) = { 0.91057460456685, -0.01210760938604, 0.0};
+Point( 273 ) = { 0.91646062035505, -0.01134951733663, 0.0};
+Point( 274 ) = { 0.92216396275101, -0.01060996200508, 0.0};
+Point( 275 ) = { 0.92768213008025, -0.00988970099293, 0.0};
+Point( 276 ) = { 0.93301270189222, -0.00918948088016, 0.0};
+Point( 277 ) = { 0.93815334002193, -0.00851003605407, 0.0};
+Point( 278 ) = { 0.94310178961561, -0.00785208747637, 0.0};
+Point( 279 ) = { 0.94785588011971, -0.00721634139577, 0.0};
+Point( 280 ) = { 0.95241352623301, -0.00660348801398, 0.0};
+Point( 281 ) = { 0.9567727288213, -0.00601420011357, 0.0};
+Point( 282 ) = { 0.96093157579425, -0.00544913165641, 0.0};
+Point( 283 ) = { 0.96488824294413, -0.00490891636171, 0.0};
+Point( 284 ) = { 0.96864099474595, -0.00439416627328, 0.0};
+Point( 285 ) = { 0.97218818511874, -0.0039054703252, 0.0};
+Point( 286 ) = { 0.97552825814758, -0.00344339291591, 0.0};
+Point( 287 ) = { 0.97865974876603, -0.00300847250009, 0.0};
+Point( 288 ) = { 0.98158128339883, -0.00260122020822, 0.0};
+Point( 289 ) = { 0.98429158056432, -0.0022221185031, 0.0};
+Point( 290 ) = { 0.98678945143658, -0.00187161988278, 0.0};
+Point( 291 ) = { 0.9890738003669, -0.0015501456389, 0.0};
+Point( 292 ) = { 0.99114362536434, -0.00125808467902, 0.0};
+Point( 293 ) = { 0.99299801853525, -0.00099579242133, 0.0};
+Point( 294 ) = { 0.99463616648149, -0.00076358976946, 0.0};
+Point( 295 ) = { 0.99605735065724, -0.00056176217462, 0.0};
+Point( 296 ) = { 0.99726094768414, -0.00039055879171, 0.0};
+Point( 297 ) = { 0.99824642962475, -0.0002501917354, 0.0};
+Point( 298 ) = { 0.99901336421414, -0.00014083544147, 0.0};
+Point( 299 ) = { 0.99956141504943, -6.262613794e-05, 0.0};
+Point( 300 ) = { 0.99989034173742, -1.566142976e-05, 0.0};
+// Point( 301 ) = { 1.0, 0.0, 0.0, msTe };
+
+Spline(1) = {Le:1}; // upper side
+Spline(2) = {1, 300:Le}; // lower side
+
+// Rotation
+If (angle != 0)
+   For i In {Te:N:1}
+       Rotate{{0, 0, 1}, {xRot, 0, 0}, -angle} {Point{i};}
+   EndFor
+EndIf
+
+// Farfield
+Point(10001) = {1+xLgt, 0, 0,msF};
+Point(10002) = {1+xLgt, yLgt, 0,msF};
+Point(10003) = {-xLgt, yLgt, 0,msF};
+Point(10004) = {-xLgt, 0, 0,msF};
+Point(10005) = {-xLgt,-yLgt, 0,msF};
+Point(10006) = {1+xLgt, -yLgt, 0,msF};
+
+Line(10001) = {10001, 10002};
+Line(10002) = {10002, 10003};
+Line(10003) = {10003, 10004};
+Line(10004) = {10004, 10005};
+Line(10005) = {10005, 10006};
+Line(10006) = {10006, 10001};
+
+// Front and wake
+Line(10007) = {Le, 10004};
+Line(10008) = {Te, 10001};
+
+// Internal field
+Line Loop(20001) = {10008, 10001, 10002, 10003, -10007, 1};
+Line Loop(20002) = {10007, 10004, 10005, 10006, -10008, 2};
+Plane Surface(30001) = {20001};
+Plane Surface(30002) = {20002};
+
+/************************* 
+ Mesh Options 
+ *************************/
+
+Mesh.Algorithm = 5; // Delaunay
+
+/************************* 
+ Physical Groups 
+ *************************/
+
+Physical Point("te") = {Te};
+Physical Line("upstream") = {10003, 10004};
+Physical Line("farfield") = {10002, 10005};
+Physical Line("downstream") = {10001};
+Physical Line("downstream") += {10006};
+Physical Line("airfoil") = {1};
+Physical Line("airfoil_") = {2};
+Physical Line("wake") = {10008};
+Physical Surface("field") = {30001};
+Physical Surface("field") += {30002};
diff --git a/blast/thesis_tests/NACA4412.geo b/blast/thesis_tests/NACA4412.geo
new file mode 100644
index 0000000000000000000000000000000000000000..b0eab3b8c543d54b8eb579262872be0bccc47fe9
--- /dev/null
+++ b/blast/thesis_tests/NACA4412.geo
@@ -0,0 +1,389 @@
+/* Airfoil NACA4412 */
+// This file was generated automatically using geoFoils v1.0 (https://github.com/Paul-Dech/geoFoils)
+// @date: 2024-09-06 10:30:47
+// @author: Paul Dechamps
+
+// Geometry
+DefineConstant[ xLgt = { 50.0, Name "Domain length (x-dir)" }  ];
+DefineConstant[ yLgt = { 50.0, Name "Domain length (y-dir)" }  ];
+
+// Mesh
+DefineConstant[ growthRatio = { 1.5, Name "Growth Ratio" }  ];
+DefineConstant[ msTe = { 0.01, Name "Airfoil TE mesh size" }  ];
+DefineConstant[ msLe = { 0.001, Name "Airfoil LE mesh size" }  ];
+
+// Rotation
+DefineConstant[ xRot = { 0.25, Name "Center of rotation" }  ];
+DefineConstant[ angle = { 0*Pi/180, Name "Angle of rotation" }  ];
+Geometry.AutoCoherence = 0; // Needed so that gmsh does not remove duplicate
+
+If (growthRatio == 1.0)
+   msF = msLe;
+Else
+   n = Log(1 - (1 - growthRatio) * xLgt / msLe) / Log(growthRatio);
+   msF = msLe * growthRatio^(n - 1);
+EndIf
+
+/**************
+ Geometry
+ **************/
+Te = 1; // trailing edge
+Le = 151; // leading edge
+
+Point( 1 ) = { 1.0, -0.0, 0.0, msTe };
+Point( 2 ) = { 0.99989241123855, 3.014386127e-05, 0.0};
+Point( 3 ) = { 0.99956968600815, 0.000120534188, 0.0};
+Point( 4 ) = { 0.9990319474757, 0.00027104730258, 0.0};
+Point( 5 ) = { 0.99827940093596, 0.00048147738708, 0.0};
+Point( 6 ) = { 0.99731233383597, 0.0007515369488, 0.0};
+Point( 7 ) = { 0.99613111580815, 0.00108085746888, 0.0};
+Point( 8 ) = { 0.99473619871106, 0.00146899023062, 0.0};
+Point( 9 ) = { 0.99312811667654, 0.00191540732356, 0.0};
+Point( 10 ) = { 0.99130748616156, 0.00241950281844, 0.0};
+Point( 11 ) = { 0.98927500600311, 0.0029805941072, 0.0};
+Point( 12 ) = { 0.98703145747381, 0.00359792340152, 0.0};
+Point( 13 ) = { 0.98457770433631, 0.00427065938282, 0.0};
+Point( 14 ) = { 0.98191469289367, 0.00499789899555, 0.0};
+Point( 15 ) = { 0.9790434520334, 0.0057786693755, 0.0};
+Point( 16 ) = { 0.97596509326219, 0.00661192990378, 0.0};
+Point( 17 ) = { 0.97268081072861, 0.00749657437701, 0.0};
+Point( 18 ) = { 0.96919188123081, 0.00843143328358, 0.0};
+Point( 19 ) = { 0.96549966420613, 0.00941527617557, 0.0};
+Point( 20 ) = { 0.96160560169982, 0.01044681412564, 0.0};
+Point( 21 ) = { 0.9575112183097, 0.01152470225783, 0.0};
+Point( 22 ) = { 0.95321812110381, 0.01264754234126, 0.0};
+Point( 23 ) = { 0.94872799950833, 0.01381388543552, 0.0};
+Point( 24 ) = { 0.94404262516275, 0.01502223457647, 0.0};
+Point( 25 ) = { 0.93916385173977, 0.01627104749146, 0.0};
+Point( 26 ) = { 0.93409361472737, 0.01755873933295, 0.0};
+Point( 27 ) = { 0.92883393117062, 0.01888368541972, 0.0};
+Point( 28 ) = { 0.92338689937118, 0.02024422397539, 0.0};
+Point( 29 ) = { 0.91775469854243, 0.02163865885395, 0.0};
+Point( 30 ) = { 0.91193958841846, 0.0230652622429, 0.0};
+Point( 31 ) = { 0.90594390881555, 0.02452227733465, 0.0};
+Point( 32 ) = { 0.89977007914474, 0.02600792095765, 0.0};
+Point( 33 ) = { 0.89342059787452, 0.02752038615941, 0.0};
+Point( 34 ) = { 0.88689804194287, 0.0290578447338, 0.0};
+Point( 35 ) = { 0.88020506611822, 0.03061844968637, 0.0};
+Point( 36 ) = { 0.87334440230893, 0.03220033763159, 0.0};
+Point( 37 ) = { 0.86631885882143, 0.03380163111704, 0.0};
+Point( 38 ) = { 0.85913131956715, 0.03542044087026, 0.0};
+Point( 39 ) = { 0.85178474321881, 0.03705486796498, 0.0};
+Point( 40 ) = { 0.84428216231661, 0.03870300590384, 0.0};
+Point( 41 ) = { 0.8366266823253, 0.04036294261635, 0.0};
+Point( 42 ) = { 0.82882148064313, 0.04203276237077, 0.0};
+Point( 43 ) = { 0.8208698055639, 0.04371054760031, 0.0};
+Point( 44 ) = { 0.81277497519349, 0.04539438064421, 0.0};
+Point( 45 ) = { 0.80454037632216, 0.04708234540561, 0.0};
+Point( 46 ) = { 0.79616946325443, 0.04877252892853, 0.0};
+Point( 47 ) = { 0.78766575659783, 0.05046302289718, 0.0};
+Point( 48 ) = { 0.77903284201233, 0.05215192506181, 0.0};
+Point( 49 ) = { 0.7702743689221, 0.05383734059537, 0.0};
+Point( 50 ) = { 0.76139404919096, 0.05551738338666, 0.0};
+Point( 51 ) = { 0.7523956557634, 0.05719017727545, 0.0};
+Point( 52 ) = { 0.74328302127238, 0.05885385723629, 0.0};
+Point( 53 ) = { 0.73406003661541, 0.06050657051748, 0.0};
+Point( 54 ) = { 0.72473064950011, 0.06214647774267, 0.0};
+Point( 55 ) = { 0.71529886296045, 0.06377175398237, 0.0};
+Point( 56 ) = { 0.70576873384447, 0.0653805898032, 0.0};
+Point( 57 ) = { 0.69614437127437, 0.06697119230266, 0.0};
+Point( 58 ) = { 0.68642993507952, 0.0685417861373, 0.0};
+Point( 59 ) = { 0.67662963420273, 0.07009061455236, 0.0};
+Point( 60 ) = { 0.66674772507997, 0.07161594042033, 0.0};
+Point( 61 ) = { 0.65678850999356, 0.07311604729632, 0.0};
+Point( 62 ) = { 0.6467563353985, 0.07458924049743, 0.0};
+Point( 63 ) = { 0.63665559022145, 0.07603384821297, 0.0};
+Point( 64 ) = { 0.62649070413184, 0.07744822265232, 0.0};
+Point( 65 ) = { 0.61626614578411, 0.0788307412361, 0.0};
+Point( 66 ) = { 0.60598642103015, 0.08017980783642, 0.0};
+Point( 67 ) = { 0.59565607110062, 0.08149385407067, 0.0};
+Point( 68 ) = { 0.585279670754, 0.08277134065319, 0.0};
+Point( 69 ) = { 0.57486182639173, 0.08401075880786, 0.0};
+Point( 70 ) = { 0.56440717413801, 0.08521063174425, 0.0};
+Point( 71 ) = { 0.55392037788261, 0.08636951619878, 0.0};
+Point( 72 ) = { 0.54340612728507, 0.08748600404161, 0.0};
+Point( 73 ) = { 0.53286913573855, 0.08855872394903, 0.0};
+Point( 74 ) = { 0.52231413829191, 0.08958634313995, 0.0};
+Point( 75 ) = { 0.51174588952817, 0.09056756917439, 0.0};
+Point( 76 ) = { 0.50116916139826, 0.09150115181065, 0.0};
+Point( 77 ) = { 0.49058874100842, 0.09238588491691, 0.0};
+Point( 78 ) = { 0.48000942836029, 0.09322060843209, 0.0};
+Point( 79 ) = { 0.46943603404272, 0.09400421036964, 0.0};
+Point( 80 ) = { 0.45887337687458, 0.09473562885721, 0.0};
+Point( 81 ) = { 0.44832628149811, 0.09541385420408, 0.0};
+Point( 82 ) = { 0.43779957592277, 0.09603793098728, 0.0};
+Point( 83 ) = { 0.42729808901973, 0.09660696014672, 0.0};
+Point( 84 ) = { 0.41682664796753, 0.09712010107875, 0.0};
+Point( 85 ) = { 0.40639007564986, 0.09757657371672, 0.0};
+Point( 86 ) = { 0.39592947995664, 0.09797339614748, 0.0};
+Point( 87 ) = { 0.38541108267205, 0.09828762422892, 0.0};
+Point( 88 ) = { 0.37494100099015, 0.09851332874847, 0.0};
+Point( 89 ) = { 0.36452442975559, 0.09865046354599, 0.0};
+Point( 90 ) = { 0.354166552819, 0.09869910762927, 0.0};
+Point( 91 ) = { 0.34387253885362, 0.0986594658543, 0.0};
+Point( 92 ) = { 0.33364753709603, 0.09853186937278, 0.0};
+Point( 93 ) = { 0.32349667302164, 0.09831677583422, 0.0};
+Point( 94 ) = { 0.31342504396654, 0.09801476933104, 0.0};
+Point( 95 ) = { 0.30343771470802, 0.09762656007535, 0.0};
+Point( 96 ) = { 0.29353971301676, 0.09715298379724, 0.0};
+Point( 97 ) = { 0.28373602519416, 0.09659500085514, 0.0};
+Point( 98 ) = { 0.27403159160908, 0.09595369504991, 0.0};
+Point( 99 ) = { 0.2644313022482, 0.09523027213524, 0.0};
+Point( 100 ) = { 0.25493999229488, 0.0944260580185, 0.0};
+Point( 101 ) = { 0.24556243775148, 0.0935424966469, 0.0};
+Point( 102 ) = { 0.23630335111985, 0.09258114757554, 0.0};
+Point( 103 ) = { 0.22716737715516, 0.09154368321516, 0.0};
+Point( 104 ) = { 0.21815908870758, 0.09043188575872, 0.0};
+Point( 105 ) = { 0.20928298266628, 0.08924764378752, 0.0};
+Point( 106 ) = { 0.20054347601977, 0.08799294855887, 0.0};
+Point( 107 ) = { 0.19194490204595, 0.08666988997911, 0.0};
+Point( 108 ) = { 0.18349150664481, 0.08528065226679, 0.0};
+Point( 109 ) = { 0.17518744482572, 0.08382750931284, 0.0};
+Point( 110 ) = { 0.16703677736067, 0.08231281974564, 0.0};
+Point( 111 ) = { 0.15904346761372, 0.08073902171046, 0.0};
+Point( 112 ) = { 0.15121137855596, 0.07910862737442, 0.0};
+Point( 113 ) = { 0.14354426997419, 0.07742421716899, 0.0};
+Point( 114 ) = { 0.13604579588039, 0.07568843378398, 0.0};
+Point( 115 ) = { 0.12871950212788, 0.07390397592773, 0.0};
+Point( 116 ) = { 0.12156882423864, 0.07207359186996, 0.0};
+Point( 117 ) = { 0.11459708544534, 0.07020007278429, 0.0};
+Point( 118 ) = { 0.10780749494979, 0.06828624590925, 0.0};
+Point( 119 ) = { 0.10120314639887, 0.06633496754691, 0.0};
+Point( 120 ) = { 0.09478701657693, 0.06434911591975, 0.0};
+Point( 121 ) = { 0.088561964313, 0.06233158390687, 0.0};
+Point( 122 ) = { 0.08253072959952, 0.06028527168161, 0.0};
+Point( 123 ) = { 0.07669593291796, 0.05821307927317, 0.0};
+Point( 124 ) = { 0.07106007476583, 0.0561178990755, 0.0};
+Point( 125 ) = { 0.06562553537807, 0.05400260832714, 0.0};
+Point( 126 ) = { 0.0603945746349, 0.05187006158595, 0.0};
+Point( 127 ) = { 0.05536933214714, 0.04972308322329, 0.0};
+Point( 128 ) = { 0.05055182750894, 0.04756445996182, 0.0};
+Point( 129 ) = { 0.04594396070723, 0.04539693348177, 0.0};
+Point( 130 ) = { 0.04154751267618, 0.04322319311989, 0.0};
+Point( 131 ) = { 0.03736414598439, 0.04104586868577, 0.0};
+Point( 132 ) = { 0.03339540564198, 0.03886752341934, 0.0};
+Point( 133 ) = { 0.0296427200142, 0.03669064711357, 0.0};
+Point( 134 ) = { 0.02610740182794, 0.03451764942567, 0.0};
+Point( 135 ) = { 0.02279064925713, 0.03235085339959, 0.0};
+Point( 136 ) = { 0.019693547073, 0.03019248922206, 0.0};
+Point( 137 ) = { 0.01681706784517, 0.02804468823378, 0.0};
+Point( 138 ) = { 0.01416207317951, 0.02590947721633, 0.0};
+Point( 139 ) = { 0.01172931497914, 0.02378877297478, 0.0};
+Point( 140 ) = { 0.00951943671502, 0.02168437723497, 0.0};
+Point( 141 ) = { 0.0075329746933, 0.01959797187312, 0.0};
+Point( 142 ) = { 0.00577035930685, 0.01753111449491, 0.0};
+Point( 143 ) = { 0.00423191625928, 0.0154852343793, 0.0};
+Point( 144 ) = { 0.00291786775042, 0.0134616288019, 0.0};
+Point( 145 ) = { 0.00182833361293, 0.01146145975074, 0.0};
+Point( 146 ) = { 0.00096333239095, 0.00948575104643, 0.0};
+Point( 147 ) = { 0.00032278235222, 0.00753538587714, 0.0};
+Point( 148 ) = { -9.349757328e-05, 0.00561110475743, 0.0};
+Point( 149 ) = { -0.00028578693461, 0.00371350391849, 0.0};
+Point( 150 ) = { -0.00025446298573, 0.00184303413598, 0.0};
+Point( 151 ) = { 0.0, 0.0, 0.0, msLe };
+Point( 152 ) = { 0.00047377951088, -0.00179917684341, 0.0};
+Point( 153 ) = { 0.00116295683575, -0.00353816611664, 0.0};
+Point( 154 ) = { 0.00206676914501, -0.00521693716817, 0.0};
+Point( 155 ) = { 0.00318435839828, -0.00683549523158, 0.0};
+Point( 156 ) = { 0.00451477224078, -0.00839388132388, 0.0};
+Point( 157 ) = { 0.00605696507259, -0.00989217225556, 0.0};
+Point( 158 ) = { 0.00780979928659, -0.01133048074951, 0.0};
+Point( 159 ) = { 0.00977204667021, -0.01270895566562, 0.0};
+Point( 160 ) = { 0.01194238996446, -0.01402778232649, 0.0};
+Point( 161 ) = { 0.01431942457289, -0.0152871829391, 0.0};
+Point( 162 ) = { 0.01690166041182, -0.01648741710627, 0.0};
+Point( 163 ) = { 0.01968752389223, -0.01762878242109, 0.0};
+Point( 164 ) = { 0.02267536002283, -0.01871161513648, 0.0};
+Point( 165 ) = { 0.02586343462276, -0.01973629090156, 0.0};
+Point( 166 ) = { 0.02924993663184, -0.02070322555574, 0.0};
+Point( 167 ) = { 0.03283298050539, -0.02161287597058, 0.0};
+Point( 168 ) = { 0.03661060868017, -0.02246574092931, 0.0};
+Point( 169 ) = { 0.04058079409755, -0.02326236203299, 0.0};
+Point( 170 ) = { 0.04474144276952, -0.024003324622, 0.0};
+Point( 171 ) = { 0.049090396373, -0.02468925870107, 0.0};
+Point( 172 ) = { 0.0536254348578, -0.02532083985588, 0.0};
+Point( 173 ) = { 0.05834427905336, -0.0258987901487, 0.0};
+Point( 174 ) = { 0.06324459325985, -0.02642387898054, 0.0};
+Point( 175 ) = { 0.06832398780899, -0.02689692390728, 0.0};
+Point( 176 ) = { 0.07358002158066, -0.02731879139673, 0.0};
+Point( 177 ) = { 0.07901020446143, -0.0276903975141, 0.0};
+Point( 178 ) = { 0.08461199973216, -0.02801270852322, 0.0};
+Point( 179 ) = { 0.09038282637194, -0.02828674139092, 0.0};
+Point( 180 ) = { 0.09632006126677, -0.02851356418254, 0.0};
+Point( 181 ) = { 0.10242104131205, -0.02869429633656, 0.0};
+Point( 182 ) = { 0.10868306539888, -0.02883010880699, 0.0};
+Point( 183 ) = { 0.11510339627529, -0.02892222406253, 0.0};
+Point( 184 ) = { 0.12167926227442, -0.02897191593207, 0.0};
+Point( 185 ) = { 0.12840785890291, -0.02898050928686, 0.0};
+Point( 186 ) = { 0.13528635028396, -0.02894937955032, 0.0};
+Point( 187 ) = { 0.14231187045071, -0.02887995202725, 0.0};
+Point( 188 ) = { 0.1494815244868, -0.02877370104505, 0.0};
+Point( 189 ) = { 0.15679238951245, -0.02863214890061, 0.0};
+Point( 190 ) = { 0.16424151551535, -0.0284568646074, 0.0};
+Point( 191 ) = { 0.17182592602743, -0.02824946243829, 0.0};
+Point( 192 ) = { 0.17954261864923, -0.02801160026098, 0.0};
+Point( 193 ) = { 0.18738856542559, -0.0277449776638, 0.0};
+Point( 194 ) = { 0.19536071307688, -0.02745133387096, 0.0};
+Point( 195 ) = { 0.20345598309167, -0.02713244544754, 0.0};
+Point( 196 ) = { 0.21167127168776, -0.02679012379581, 0.0};
+Point( 197 ) = { 0.22000344964929, -0.02642621244568, 0.0};
+Point( 198 ) = { 0.22844936204907, -0.0260425841433, 0.0};
+Point( 199 ) = { 0.23700582786584, -0.02564113774339, 0.0};
+Point( 200 ) = { 0.24566963950702, -0.02522379491189, 0.0};
+Point( 201 ) = { 0.25443756224852, -0.0247924966469, 0.0};
+Point( 202 ) = { 0.2633063336034, -0.02434919962723, 0.0};
+Point( 203 ) = { 0.27227266263194, -0.02389587239897, 0.0};
+Point( 204 ) = { 0.28133322920599, -0.02343449141176, 0.0};
+Point( 205 ) = { 0.29048468324076, -0.0229670369176, 0.0};
+Point( 206 ) = { 0.29972364390744, -0.02249548874602, 0.0};
+Point( 207 ) = { 0.30904669883987, -0.02202182197065, 0.0};
+Point( 208 ) = { 0.31845040334878, -0.02154800248296, 0.0};
+Point( 209 ) = { 0.32793127965654, -0.02107598248989, 0.0};
+Point( 210 ) = { 0.33748581616539, -0.02060769595298, 0.0};
+Point( 211 ) = { 0.34711046677144, -0.02014505398712, 0.0};
+Point( 212 ) = { 0.35680165023653, -0.01968994023766, 0.0};
+Point( 213 ) = { 0.36655574962914, -0.0192442062552, 0.0};
+Point( 214 ) = { 0.37636911184499, -0.01880966688749, 0.0};
+Point( 215 ) = { 0.3862380472173, -0.01838809570843, 0.0};
+Point( 216 ) = { 0.3961588292256, -0.01798122050393, 0.0};
+Point( 217 ) = { 0.40622860976441, -0.0175854198957, 0.0};
+Point( 218 ) = { 0.41640460531637, -0.01718145197847, 0.0};
+Point( 219 ) = { 0.42661888241785, -0.01676846235833, 0.0};
+Point( 220 ) = { 0.43686719051292, -0.01634766021005, 0.0};
+Point( 221 ) = { 0.44714525523424, -0.01592023277775, 0.0};
+Point( 222 ) = { 0.4574487797931, -0.01548734242009, 0.0};
+Point( 223 ) = { 0.46777344642796, -0.01505012378803, 0.0};
+Point( 224 ) = { 0.47811491791051, -0.01460968114784, 0.0};
+Point( 225 ) = { 0.48846883910822, -0.0141670858612, 0.0};
+Point( 226 ) = { 0.49883083860174, -0.01372337403287, 0.0};
+Point( 227 ) = { 0.50919653035519, -0.01327954433572, 0.0};
+Point( 228 ) = { 0.51956151543729, -0.01283655602144, 0.0};
+Point( 229 ) = { 0.52992138379076, -0.01239532712428, 0.0};
+Point( 230 ) = { 0.54027171604725, -0.0119567328637, 0.0};
+Point( 231 ) = { 0.55060808538504, -0.01152160425101, 0.0};
+Point( 232 ) = { 0.5609260594263, -0.01109072690322, 0.0};
+Point( 233 ) = { 0.57122120217068, -0.01066484006669, 0.0};
+Point( 234 ) = { 0.5814890759621, -0.01024463585139, 0.0};
+Point( 235 ) = { 0.5917252434851, -0.00983075867568, 0.0};
+Point( 236 ) = { 0.60192526978761, -0.00942380492007, 0.0};
+Point( 237 ) = { 0.61208472432654, -0.00902432278747, 0.0};
+Point( 238 ) = { 0.62219918303302, -0.00863281236587, 0.0};
+Point( 239 ) = { 0.63226423039382, -0.00824972588885, 0.0};
+Point( 240 ) = { 0.64227546154597, -0.00787546818758, 0.0};
+Point( 241 ) = { 0.65222848438139, -0.00751039732757, 0.0};
+Point( 242 ) = { 0.66211892165862, -0.00715482542213, 0.0};
+Point( 243 ) = { 0.67194241311909, -0.00680901961363, 0.0};
+Point( 244 ) = { 0.68169461760516, -0.00647320321303, 0.0};
+Point( 245 ) = { 0.69137121517773, -0.00614755698733, 0.0};
+Point( 246 ) = { 0.70096790923133, -0.00583222058381, 0.0};
+Point( 247 ) = { 0.71048042860462, -0.00552729407969, 0.0};
+Point( 248 ) = { 0.71990452968482, -0.00523283964504, 0.0};
+Point( 249 ) = { 0.72923599850446, -0.00494888330668, 0.0};
+Point( 250 ) = { 0.73847065282933, -0.00467541680024, 0.0};
+Point( 251 ) = { 0.7476043442366, -0.00441239949767, 0.0};
+Point( 252 ) = { 0.75663296018217, -0.00415976039717, 0.0};
+Point( 253 ) = { 0.7655524260569, -0.00391740016254, 0.0};
+Point( 254 ) = { 0.77435870723101, -0.0036851931991, 0.0};
+Point( 255 ) = { 0.78304781108661, -0.00346298975345, 0.0};
+Point( 256 ) = { 0.79161578903804, -0.00325061802461, 0.0};
+Point( 257 ) = { 0.80005873854021, -0.00304788627435, 0.0};
+Point( 258 ) = { 0.80837280508482, -0.00285458492509, 0.0};
+Point( 259 ) = { 0.81655418418479, -0.00267048863401, 0.0};
+Point( 260 ) = { 0.82459912334698, -0.00249535833271, 0.0};
+Point( 261 ) = { 0.83250392403356, -0.00232894322245, 0.0};
+Point( 262 ) = { 0.84026494361207, -0.00217098271533, 0.0};
+Point( 263 ) = { 0.84787859729455, -0.00202120831296, 0.0};
+Point( 264 ) = { 0.85534136006565, -0.00187934541442, 0.0};
+Point( 265 ) = { 0.86264976859998, -0.00174511504666, 0.0};
+Point( 266 ) = { 0.86980042316846, -0.00161823551075, 0.0};
+Point( 267 ) = { 0.87678998953353, -0.0014984239387, 0.0};
+Point( 268 ) = { 0.88361520083292, -0.0013853977562, 0.0};
+Point( 269 ) = { 0.89027285945132, -0.0012788760475, 0.0};
+Point( 270 ) = { 0.89675983887946, -0.00117858081973, 0.0};
+Point( 271 ) = { 0.9030730855594, -0.0010842381645, 0.0};
+Point( 272 ) = { 0.90920962071525, -0.00099557931592, 0.0};
+Point( 273 ) = { 0.91516654216767, -0.00091234160451, 0.0};
+Point( 274 ) = { 0.92094102613083, -0.0008342693078, 0.0};
+Point( 275 ) = { 0.92653032898989, -0.00076111439884, 0.0};
+Point( 276 ) = { 0.93193178905707, -0.00069263719483, 0.0};
+Point( 277 ) = { 0.93714282830409, -0.00062860690852, 0.0};
+Point( 278 ) = { 0.94216095406846, -0.00056880210617, 0.0};
+Point( 279 ) = { 0.94698376073108, -0.00051301107591, 0.0};
+Point( 280 ) = { 0.95160893136221, -0.0004610321113, 0.0};
+Point( 281 ) = { 0.9560342393329, -0.00041267371541, 0.0};
+Point( 282 ) = { 0.96025754988868, -0.00036775473076, 0.0};
+Point( 283 ) = { 0.96427682168212, -0.00032610440146, 0.0};
+Point( 284 ) = { 0.96809010826108, -0.00028756237373, 0.0};
+Point( 285 ) = { 0.97169555950887, -0.00025197864134, 0.0};
+Point( 286 ) = { 0.97509142303296, -0.00021921344297, 0.0};
+Point( 287 ) = { 0.97827604549867, -0.00018913711816, 0.0};
+Point( 288 ) = { 0.98124787390399, -0.00016162992884, 0.0};
+Point( 289 ) = { 0.98400545679232, -0.00013658185356, 0.0};
+Point( 290 ) = { 0.98654744539935, -0.00011389236091, 0.0};
+Point( 291 ) = { 0.9888725947307, -9.347016913e-05, 0.0};
+Point( 292 ) = { 0.99097976456712, -7.52329982e-05, 0.0};
+Point( 293 ) = { 0.99286792039397, -5.910732061e-05, 0.0};
+Point( 294 ) = { 0.99453613425193, -4.50281168e-05, 0.0};
+Point( 295 ) = { 0.99598358550633, -3.293864056e-05, 0.0};
+Point( 296 ) = { 0.99720956153231, -2.279019959e-05, 0.0};
+Point( 297 ) = { 0.99821345831355, -1.454195569e-05, 0.0};
+Point( 298 ) = { 0.99899478095257, -8.16074861e-06, 0.0};
+Point( 299 ) = { 0.99955314409071, -3.62094713e-06, 0.0};
+Point( 300 ) = { 0.9998882722363, -9.0433012e-07, 0.0};
+// Point( 301 ) = { 1.0, 0.0, 0.0, msTe };
+
+Spline(1) = {Le:1}; // upper side
+Spline(2) = {1, 300:Le}; // lower side
+
+// Rotation
+If (angle != 0)
+   For i In {Te:N:1}
+       Rotate{{0, 0, 1}, {xRot, 0, 0}, -angle} {Point{i};}
+   EndFor
+EndIf
+
+// Farfield
+Point(10001) = {1+xLgt, 0, 0,msF};
+Point(10002) = {1+xLgt, yLgt, 0,msF};
+Point(10003) = {-xLgt, yLgt, 0,msF};
+Point(10004) = {-xLgt, 0, 0,msF};
+Point(10005) = {-xLgt,-yLgt, 0,msF};
+Point(10006) = {1+xLgt, -yLgt, 0,msF};
+
+Line(10001) = {10001, 10002};
+Line(10002) = {10002, 10003};
+Line(10003) = {10003, 10004};
+Line(10004) = {10004, 10005};
+Line(10005) = {10005, 10006};
+Line(10006) = {10006, 10001};
+
+// Front and wake
+Line(10007) = {Le, 10004};
+Line(10008) = {Te, 10001};
+
+// Internal field
+Line Loop(20001) = {10008, 10001, 10002, 10003, -10007, 1};
+Line Loop(20002) = {10007, 10004, 10005, 10006, -10008, 2};
+Plane Surface(30001) = {20001};
+Plane Surface(30002) = {20002};
+
+/************************* 
+ Mesh Options 
+ *************************/
+
+Mesh.Algorithm = 5; // Delaunay
+
+/************************* 
+ Physical Groups 
+ *************************/
+
+Physical Point("te") = {Te};
+Physical Line("upstream") = {10003, 10004};
+Physical Line("farfield") = {10002, 10005};
+Physical Line("downstream") = {10001};
+Physical Line("downstream") += {10006};
+Physical Line("airfoil") = {1};
+Physical Line("airfoil_") = {2};
+Physical Line("wake") = {10008};
+Physical Surface("field") = {30001};
+Physical Surface("field") += {30002};
diff --git a/blast/thesis_tests/naca0012_2D.py b/blast/thesis_tests/naca0012_2D.py
new file mode 100644
index 0000000000000000000000000000000000000000..7b3cbaf590b22718ab72426b3d98e8eef5e1c765
--- /dev/null
+++ b/blast/thesis_tests/naca0012_2D.py
@@ -0,0 +1,173 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Imports.
+
+import blast.utils as viscUtils
+import numpy as np
+import matplotlib.pyplot as plt
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def cfgInviscid(nthrds, verb):
+    import os.path
+    # Parameters
+    return {
+    # Options
+    'scenario' : 'aerodynamic',
+    'task' : 'analysis',
+    'Threads' : nthrds, # number of threads
+    'Verb' : verb, # verbosity
+    # Model (geometry or mesh)
+    'File' : os.path.dirname(os.path.abspath(__file__)) + '/NACA0012.geo', # Input file containing the model
+    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msTe' : 0.01, 'msLe' : 0.001, 'growthRatio' : 1.12}, # parameters for input file model
+    'Dim' : 2, # problem dimension
+    'Format' : 'gmsh', # save format (vtk or gmsh)
+    # Markers
+    'Fluid' : 'field', # name of physical group containing the fluid
+    'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+    'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
+    'Wake' : 'wake', # LIST of names of physical group containing the wake
+    'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
+    'Te' : 'te', # LIST of names of physical group containing the trailing edge
+    'dbc' : True,
+    'Upstream' : 'upstream',
+    # Freestream
+    'M_inf' : .4, # freestream Mach number
+    'AoA' : 5., # freestream angle of attack
+    # Geometry
+    'S_ref' : 1., # reference surface length
+    'c_ref' : 1., # reference chord length
+    'x_ref' : .25, # reference point for moment computation (x)
+    'y_ref' : 0.0, # reference point for moment computation (y)
+    'z_ref' : 0.0, # reference point for moment computation (z)
+    # Numerical
+    'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
+    'G_fill' : 2, # fill-in factor for GMRES preconditioner
+    'G_tol' : 1e-5, # tolerance for GMRES
+    'G_restart' : 50, # restart for GMRES
+    'Rel_tol' : 1e-14, # relative tolerance on solver residual
+    'Abs_tol' : 1e-12, # absolute tolerance on solver residual
+    'Max_it' : 50 # solver maximum number of iterations
+    }
+
+def cfgBlast(verb):
+    return {
+        'Re' : 50e3,                 # Freestream Reynolds number
+        'Minf' : .4,               # Freestream Mach number (used for the computation of the time step only)
+        'CFL0' : 1,                 # Inital CFL number of the calculation
+        'Verb': verb,               # Verbosity level of the solver
+        'couplIter': 500,            # Maximum number of iterations
+        'couplTol' : 1e-4,          # Tolerance of the VII methodology
+        'iterPrint': 1,             # int, number of iterations between outputs
+        'resetInv' : True,          # bool, flag to reset the inviscid calculation at every iteration.
+        'sections' : [0],           # List of sections for boundary layer calculation
+        'xtrF' : [-1, -1],      # Forced transition location
+        'interpolator' : 'Matching' # Interpolator for the coupling
+    }
+
+def main():
+    # Timer.
+
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    args = parseargs()
+    icfg = cfgInviscid(args.k, args.verb)
+    vcfg = cfgBlast(args.verb)
+
+    tms['pre'].start()
+    coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
+    tms['pre'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    aeroCoeffs = coupler.run()
+    tms['solver'].stop()
+
+    # Display results.
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
+
+     # Write results to file.
+    vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all')
+    tms['total'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    print('SOLVERS statistics')
+    print(coupler.tms)
+    
+
+    
+    # Show results
+    if not args.nogui:
+        iCp = viscUtils.read('Cp_inviscid.dat')
+        vCp = viscUtils.read('Cp_viscous.dat')
+        #xfoilCp = viscUtils.read('../../blast/thesis_tests/cp_naca4412f.dat', delim=None, skip = 1) #
+        #xfoilCpInv = viscUtils.read('../../blast/thesis_tests/cp_naca4412invf.dat', delim=None, skip = 1) #
+        plotcp = {'curves': [np.column_stack((vCp[:,0], vCp[:,3])),
+                             np.column_stack((iCp[:,0], iCp[:,3]))], #np.column_stack((xfoilCp[:,0], xfoilCp[:,1])), np.column_stack((xfoilCpInv[:,0], xfoilCpInv[:,1]))],
+                    'labels': ['Blaster (VII)', 'DART (inviscid)'],#'XFoil (viscous)','XFoil (inviscid)'],
+                    'lw': [3, 2, 3, 2],
+                    'color': ['firebrick','darkblue','firebrick', 'darkblue'],
+                    'ls': ['-', '-','--','--'],
+                    'reverse': True,
+                    'xlim':[0, 1],
+                    'yreverse': True,
+                    'legend': True,
+                    'xlabel': '$x/c$',
+                    'ylabel': '$c_p$'
+                    }
+        viscUtils.plot(plotcp)
+  
+"""
+# Initialize an empty list to store x-coordinates
+    x_coordinate = []
+    y_coordinate = []
+   #mach = []
+
+    with open("boundary_coordinates.txt", "r") as file:
+        for line in file:
+            coordinates = line.strip().split()  # Split the line into x and y values
+            if coordinates:  # Ensure the line is not empty
+                x_coordinate.append(coordinates[0])  # Append only the x-coordinate
+                y_coordinate.append(coordinates[1])
+            else:
+                break  # Stop reading if an empty line is encountered (optional)
+
+    l = np.linspace(0, 100,np.size(x_coordinate), endpoint=True)
+
+with open("viscous_input_mach.txt", "r") as file:
+    for line in file:
+        m = line.strip().split() 
+        if m:  # Ensure the line is not empty
+            mach.append(m[0]) 
+            break  # Stop reading if an empty line is encountered (optional)
+  
+
+
+    fig, ax = plt.subplots()
+    ax.plot(x_coordinate,l)  
+    #plt.title('Mach at boundary')
+    #plt.xlabel(r'x/c [-]')
+    #plt.ylabel(r'M [-]')
+    #ax.set_xlim([0,1])
+    #ax.set_ylim([0,0.12])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    plt.show()  # Display the plot
+"""
+
+# eof
+print('')
+
+if __name__ == "__main__":
+    main()
diff --git a/blast/thesis_tests/naca4412_2D.py b/blast/thesis_tests/naca4412_2D.py
new file mode 100644
index 0000000000000000000000000000000000000000..a34492113cc793bc834e97cc4f13e5ce151ce1ce
--- /dev/null
+++ b/blast/thesis_tests/naca4412_2D.py
@@ -0,0 +1,168 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Imports.
+
+import blast.utils as viscUtils
+import numpy as np
+import matplotlib.pyplot as plt
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def cfgInviscid(nthrds, verb):
+    import os.path
+    # Parameters
+    return {
+    # Options
+    'scenario' : 'aerodynamic',
+    'task' : 'analysis',
+    'Threads' : nthrds, # number of threads
+    'Verb' : verb, # verbosity
+    # Model (geometry or mesh)
+    'File' : os.path.dirname(os.path.abspath(__file__)) + '/NACA4412.geo', # Input file containing the model
+    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msTe' : 0.01, 'msLe' : 0.001, 'growthRatio' : 1.2}, # parameters for input file model
+    'Dim' : 2, # problem dimension
+    'Format' : 'gmsh', # save format (vtk or gmsh)
+    # Markers
+    'Fluid' : 'field', # name of physical group containing the fluid
+    'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+    'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
+    'Wake' : 'wake', # LIST of names of physical group containing the wake
+    'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
+    'Te' : 'te', # LIST of names of physical group containing the trailing edge
+    'dbc' : True,
+    'Upstream' : 'upstream',
+    # Freestream
+    'M_inf' : .3, # freestream Mach number
+    'AoA' : 0., # freestream angle of attack
+    # Geometry
+    'S_ref' : 1., # reference surface length
+    'c_ref' : 1., # reference chord length
+    'x_ref' : .25, # reference point for moment computation (x)
+    'y_ref' : 0.0, # reference point for moment computation (y)
+    'z_ref' : 0.0, # reference point for moment computation (z)
+    # Numerical
+    'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
+    'G_fill' : 2, # fill-in factor for GMRES preconditioner
+    'G_tol' : 1e-5, # tolerance for GMRES
+    'G_restart' : 50, # restart for GMRES
+    'Rel_tol' : 1e-6, # relative tolerance on solver residual
+    'Abs_tol' : 1e-8, # absolute tolerance on solver residual
+    'Max_it' : 20 # solver maximum number of iterations
+    }
+
+def cfgBlast(verb):
+    return {
+        'Re' : 1e7,                 # Freestream Reynolds number
+        'Minf' : .3,               # Freestream Mach number (used for the computation of the time step only)
+        'CFL0' : 1,                 # Inital CFL number of the calculation
+        'Verb': verb,               # Verbosity level of the solver
+        'couplIter': 100,            # Maximum number of iterations
+        'couplTol' : 1e-4,          # Tolerance of the VII methodology
+        'iterPrint': 5,             # int, number of iterations between outputs
+        'resetInv' : True,          # bool, flag to reset the inviscid calculation at every iteration.
+        'sections' : [0],           # List of sections for boundary layer calculation
+        'xtrF' : [0, 0],       # Forced transition location
+        'interpolator' : 'Matching' # Interpolator for the coupling
+    }
+
+def main():
+    # Timer.
+
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    args = parseargs()
+    icfg = cfgInviscid(args.k, args.verb)
+    vcfg = cfgBlast(args.verb)
+
+    tms['pre'].start()
+    coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
+    tms['pre'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    aeroCoeffs = coupler.run()
+    tms['solver'].stop()
+
+    # Display results.
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
+
+     # Write results to file.
+    vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all')
+    tms['total'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    print('SOLVERS statistics')
+    print(coupler.tms)
+    
+
+
+    # Show results
+    if not args.nogui:
+        iCp = viscUtils.read('Cp_inviscid.dat')
+        vCp = viscUtils.read('Cp_viscous.dat')
+        xfoilCp = viscUtils.read('../../blast/thesis_tests/cp_naca4412f.dat', delim=None, skip = 1) #
+        xfoilCpInv = viscUtils.read('../../blast/thesis_tests/cp_naca4412invf.dat', delim=None, skip = 1) #
+        plotcp = {'curves': [np.column_stack((vCp[:,0], vCp[:,3])),
+                             np.column_stack((xfoilCp[:,0], xfoilCp[:,1])), 
+                             np.column_stack((iCp[:,0], iCp[:,3])),
+                             np.column_stack((xfoilCpInv[:,0], xfoilCpInv[:,1]))],
+                    'labels': ['Blaster (VII)', 'XFoil (viscous)','DART (inviscid)','XFoil (inviscid)'],
+                    'lw': [3, 2, 3, 2],
+                    'color': ['firebrick','darkblue','firebrick', 'darkblue'],
+                    'ls': ['-', '-','--','--'],
+                    'reverse': True,
+                    'xlim':[0, 1],
+                    'yreverse': True,
+                    'legend': True,
+                    'xlabel': '$x/c$',
+                    'ylabel': '$c_p$'
+                    }
+        viscUtils.plot(plotcp)
+
+
+# Initialize an empty list to store x-coordinates
+    x_coordinate = []
+    mach = []
+"""
+    with open("boundary_coordinates.txt", "r") as file:
+        for line in file:
+            x = line.strip().split()  # Split the line into x and y values
+            if x:  # Ensure the line is not empty
+                x_coordinate.append(x[0])  # Append only the x-coordinate
+            else:
+                break  # Stop reading if an empty line is encountered (optional)
+
+    with open("viscous_input_mach.txt", "r") as file:
+        for line in file:
+            m = line.strip().split() 
+            if m:  # Ensure the line is not empty
+                mach.append(m[0]) 
+                break  # Stop reading if an empty line is encountered (optional)
+
+
+    fig, ax = plt.subplots()
+    ax.plot(mach,x_coordinate[458:499])  
+    plt.title('Mach at boundary')
+    plt.xlabel(r'x/c [-]')
+    plt.ylabel(r'M [-]')
+    ax.set_xticks([])
+    ax.set_yticks([])
+    plt.show()  # Display the plot
+"""
+
+# eof
+print('')
+
+if __name__ == "__main__":
+    main()