diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index e0d5813bfb4f6b51c38e086d2567b7cf5db612e1..887e692cf6f87f5435a20da24803f444b6da4683 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -89,19 +89,19 @@ class Coupler:
       error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
       
       if error  <= self.couplTol and exitCode==0:
-        print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
-        print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}\n'.format(couplIter,
+        print('{:>5s}| {:>10s} {:>10s} | {:>10s} {:>8s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
+        print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.2f}\n'.format(couplIter,
         self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)),ccolors.ANSI_RESET)
         return 0
       
       cdPrev = self.vSolver.Cdt
 
 
-      if couplIter%5==0:
+      if couplIter%1==0:
         xtr = [1, 1]
         xtr[0]=self.vSolver.Getxtr(0,0)
         xtr[1]=self.vSolver.Getxtr(0,1)
-        print(' {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}'.format(couplIter,
+        print(' {:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.3f}'.format(couplIter,
         self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)))
 
       couplIter += 1
@@ -157,6 +157,9 @@ class Coupler:
 
     self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()]
     self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()]
+    """
+    plt.plot(self.group[0].u)
+    plt.show()"""
 
     for n in range(0, len(self.group)):
       for i in range(0, self.group[n].nE):
@@ -209,6 +212,31 @@ class Coupler:
       self.vSolver.SetDeltaStarExt(0, 1, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])    
       self.vSolver.SetDeltaStarExt(0, 2, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
 
+    """plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['vx'][:fMeshDict['bNodes'][1]], 'x')
+    plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vx'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], 'x')
+    plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][2]:], fMeshDict['vx'][fMeshDict['bNodes'][2]:], 'x')
+    plt.title('Ux = f(x)')
+    plt.show()
+    plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['vy'][:fMeshDict['bNodes'][1]], 'x')
+    plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vy'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], 'x')
+    plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][2]:], fMeshDict['vy'][fMeshDict['bNodes'][2]:], 'x')
+    plt.title('Uy = f(x)')
+    plt.show()
+    plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['vz'][:fMeshDict['bNodes'][1]], 'x')
+    plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vz'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], 'x')
+    plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][2]:], fMeshDict['vz'][fMeshDict['bNodes'][2]:], 'x')
+    plt.title('Uz = f(x)')
+    plt.show()"""
+
+    """plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['vy'][:fMeshDict['bNodes'][1]], 'x-')
+    plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vy'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], 'x-')
+    #plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][2]:], fMeshDict['vy'][fMeshDict['bNodes'][2]:], 'x-')
+    plt.xlim([0.4, 1.2])
+    plt.ylim([-0.3, 0.2])
+    plt.title('Uy = f(x)')
+    plt.show()"""
+
+
     self.nElmUpperPrev = len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]])
     self.vSolver.SetVelocities(0, 0, fMeshDict['vx'][:fMeshDict['bNodes'][1]], fMeshDict['vy'][:fMeshDict['bNodes'][1]], fMeshDict['vz'][:fMeshDict['bNodes'][1]])
     self.vSolver.SetVelocities(0, 1, fMeshDict['vx'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vy'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vz'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
diff --git a/dart/pyVII/vCoupler2.py b/dart/pyVII/vCoupler2.py
new file mode 100644
index 0000000000000000000000000000000000000000..d0f24b4efa984dc0dabbcb4789f8f8d0bfc44925
--- /dev/null
+++ b/dart/pyVII/vCoupler2.py
@@ -0,0 +1,130 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## VII Coupler
+# Paul Dechamps
+
+from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
+
+from fwk.coloring import ccolors
+import math
+
+import fwk
+import dart.pyVII.vUtils as viscU
+import numpy as np
+from matplotlib import pyplot as plt
+
+class Coupler:
+  def __init__(self, _iSolverAPI, vSolver) -> None:
+    self.iSolverAPI = _iSolverAPI
+    self.vSolver = vSolver
+
+    self.maxCouplIter = 150
+    self.couplTol = 1e-4
+
+    self.tms=fwk.Timers()
+    self.resetInv = True
+    pass
+
+  def Run(self):
+    self.tms['tot'].start()
+
+    # Initilize meshes in c++ objects
+
+    dragIter = []
+
+    self.nElmUpperPrev = 0
+    couplIter = 0
+    cdPrev = 0
+    print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'
+    .format(self.iSolverAPI.GetAlpha()*180/math.pi, self.iSolverAPI.GetMinf(), self.vSolver.GetRe()/1e6))
+    print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
+    print(' --------------------------------------------------------')
+    while couplIter < self.maxCouplIter:
+
+      # Run inviscid solver
+      self.tms['inviscid'].start()
+      if self.resetInv:
+        self.iSolverAPI.ResetSolver()
+      self.iSolverAPI.RunSolver()
+      self.tms['inviscid'].stop()
+
+      # Extract inviscid boundary
+      self.tms['processing'].start()
+      self.iSolverAPI.GetInvBC(self.vSolver, couplIter)
+      self.tms['processing'].stop()
+
+      # Run viscous solver
+      self.tms['viscous'].start()
+      exitCode = self.vSolver.Run(couplIter)
+      self.tms['viscous'].stop()
+      # print('Viscous ops terminated with exit code ', exitCode)
+
+      self.tms['processing'].start()
+      self.iSolverAPI.SetBlowingVel(self.vSolver)
+      self.tms['processing'].stop()
+
+      dragIter.append(self.vSolver.Cdt)
+      error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
+      
+      if error  <= self.couplTol:
+
+        print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.3f} {:>6.3f} {:>8.2f}\n'.format(couplIter,
+        self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)),ccolors.ANSI_RESET)
+        self.tms['tot'].stop()
+        
+        return dragIter
+      
+      cdPrev = self.vSolver.Cdt
+
+
+      print('- {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.3f} {:>6.3f} {:>8.2f}'.format(couplIter,
+      self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)))
+
+      couplIter += 1
+      """if couplIter == 5:
+        quit()"""
+      self.tms['processing'].stop()
+    else:
+      print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
+      self.tms['tot'].stop()
+
+      return dragIter
+
+
+  def RunPolar(self, alphaSeq):
+    
+    import math
+
+    self.iSolver.printOn = False
+    self.vSolver.printOn = False
+
+    coeffTable = np.empty((3,len(alphaSeq)))
+    coeffTable[0,:] = alphaSeq
+
+    for ind, alpha in enumerate(alphaSeq):
+      self.iSolver.pbl.update(alpha*math.pi/180)
+
+      self.Run()
+
+      # Collect Coeff
+
+      coeffTable[1,ind] = self.iSolver.Cl
+      coeffTable[2,ind] = self.vSolver.Cdt
+
+    return coeffTable
diff --git a/dart/pyVII/vInterpolator.py b/dart/pyVII/vInterpolator.py
index b259a98e00d68d58d69a21b27ced9acd6ed9da36..5542ce971bfd5e4316ad26cf2aa1477ffee37250 100644
--- a/dart/pyVII/vInterpolator.py
+++ b/dart/pyVII/vInterpolator.py
@@ -1,26 +1,41 @@
-import fwk
-
 import numpy as np
 from matplotlib import pyplot as plt
+from fwk.coloring import ccolors
+
+# 3D spline with Bézier curves
 from scipy.interpolate import bisplrep
 from scipy.interpolate import bisplev
+
+# 2D spline with Bézier curves
+from scipy.interpolate import splprep
+from scipy.interpolate import splev
+from scipy.integrate import cumtrapz
+
+from scipy.interpolate import NearestNDInterpolator
+from scipy.interpolate import griddata
+from scipy.interpolate import CloughTocher2DInterpolator
+
+from dart.pyVII import vRbf as RbfInterpol
+
 from scipy.interpolate import interp1d
+
 from scipy.interpolate import RBFInterpolator
 
-class Interpolator:
+import fwk
+
+
 
-    def __init__(self, _dartSolver, bnd_wing, bnd_wake, nSections, nDim, nPoints=[500, 25], eps=0.05, k=[5, 5]) -> None:
+class Interpolator:
+    def __init__(self, _dartSolver, bnd_wing, bnd_wake, _cfg):
         self.iSolver = _dartSolver
 
         self.bndWing = bnd_wing
         self.bndWake = bnd_wake
+        self.invStruct, self.viscStruct = self.__GetWing(bnd_wing, bnd_wake, _cfg)
+        self.config = _cfg
 
-        self.nDim = nDim
-
-        self.Sections, self.Wing, self.Wake = self.__GetWing(nSections, nDim, nPoints, eps, k)
-        self.xxPrev = [[np.zeros(len(self.Sections[iSec][iReg][:,0])) for iReg in range(len(self.Sections[iSec]))] for iSec in range(len(self.Sections))]
-        self.DeltaStarPrev = [[np.zeros(len(self.Sections[iSec][iReg][:,0])) for iReg in range(len(self.Sections[iSec]))] for iSec in range(len(self.Sections))]
-        self.interpMode = 2 # 0 For Rbf, 1 for 1d interpolation.
+        self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
+        self.DeltaStarPrev = [[np.zeros(len(self.viscStruct[iSec][iReg][:,0])) for iReg in range(2)] for iSec in range(len(self.viscStruct))]
         self.tms = fwk.Timers()
 
     def GetInvBC(self, vSolver, couplIter):
@@ -51,153 +66,178 @@ class Interpolator:
             U_wk[iNode,2] = self.iSolver.U[n.row][2]
             M_wk[iNode] = self.iSolver.M[n.row]
             Rho_wk[iNode] = self.iSolver.rho[n.row]
-        
-        self.tms['interpolation'].start()
-        if self.interpMode == 0:
-            Ux = self.__RealToSections(U_wg[:,0], U_wk[:,0])
-            Uy = self.__RealToSections(U_wg[:,1], U_wk[:,1])
-            Uz = self.__RealToSections(U_wg[:,2], U_wk[:,2])
-
-            Me = self.__RealToSections(M_wg, M_wk)
-            Rho = self.__RealToSections(Rho_wg, Rho_wk)
-        elif self.interpMode == 1:
-            Ux = self.__RbfToSections(U_wg[:,0], U_wk[:,0])
-            Uy = self.__RbfToSections(U_wg[:,1], U_wk[:,1])
-            Uz = self.__RbfToSections(U_wg[:,2], U_wk[:,2])
-
-            Me = self.__RbfToSections(M_wg, M_wk)
-            Rho = self.__RbfToSections(Rho_wg, Rho_wk)
-        elif self.interpMode == 2:
-            Ux = self.__InterpToSections(U_wg[:,0], U_wk[:,0])
-            Uy = self.__InterpToSections(U_wg[:,1], U_wk[:,1])
-            Uz = self.__InterpToSections(U_wg[:,2], U_wk[:,2])
-
-            Me = self.__InterpToSections(M_wg, M_wk)
-            Rho = self.__InterpToSections(Rho_wg, Rho_wk)
-        else:
-            raise RuntimeError('Unknown interpolation mode ', self.interpMode)
-        self.tms['interpolation'].stop()
+
+        if self.config['nDim'] == 3:  # x and y inverted in 3d
+            U_wg[:, [1,2]] = U_wg[:,[2,1]]
+            U_wk[:, [1,2]] = U_wk[:,[2,1]]
+
+        Ux = self.__RbfToSections(U_wg[:,0], U_wk[:,0])
+        Uy = self.__RbfToSections(U_wg[:,1], U_wk[:,1])
+        Uz = self.__RbfToSections(U_wg[:,2], U_wk[:,2])
+
+        Me = self.__RbfToSections(M_wg, M_wk)
+        Rho = self.__RbfToSections(Rho_wg, Rho_wk)
 
         # Find stagnation point
-        for iSec in range(len(self.Sections)):
+        for iSec in range(len(self.viscStruct)):
 
-            idxStgUp = np.argmin(np.sqrt(Ux[iSec][0]**2+Uy[iSec][0]**2))
-            idxStgLw = np.argmin(np.sqrt(Ux[iSec][1]**2+Uy[iSec][1]**2))
+            #if couplIter == 0:
+            self.stgPt = np.argmin(Ux[iSec][0]**2+Uy[iSec][0]**2)
 
-            if np.sqrt(Ux[iSec][0][idxStgUp]**2+Uy[iSec][0][idxStgUp]**2) <= np.sqrt(Ux[iSec][1][idxStgLw]**2+Uy[iSec][1][idxStgLw]**2):
-                reg = 0
-                idx = idxStgUp
-            else:
-                reg = 1
-                idx = idxStgLw
-
-            xUp, xLw = self.__Remesh(self.Sections[iSec][0][:,0], self.Sections[iSec][1][:,0], idx, reg)
-            yUp, yLw = self.__Remesh(self.Sections[iSec][0][:,1], self.Sections[iSec][1][:,1], idx, reg)
-            zUp, zLw = self.__Remesh(self.Sections[iSec][0][:,2], self.Sections[iSec][1][:,2], idx, reg)
-
-            vSolver.SetGlobMesh(iSec, 0, xUp, zUp, yUp)
-            vSolver.SetGlobMesh(iSec, 1, xLw, zLw, yLw)
-            vSolver.SetGlobMesh(iSec, 2, self.Sections[iSec][2][:,0],
-                                         self.Sections[iSec][2][:,1],
-                                         self.Sections[iSec][2][:,2])
+            xUp, xLw = self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt)
+            yUp, yLw = self.__Remesh(self.viscStruct[iSec][0][:,1], self.stgPt)
+            zUp, zLw = self.__Remesh(self.viscStruct[iSec][0][:,2], self.stgPt)
+
+            vSolver.SetGlobMesh(iSec, 0, xUp, yUp, zUp)
+            vSolver.SetGlobMesh(iSec, 1, xLw, yLw, zLw)
+            vSolver.SetGlobMesh(iSec, 2, self.viscStruct[iSec][1][:,0],
+                                         self.viscStruct[iSec][1][:,1],
+                                         self.viscStruct[iSec][1][:,2])
+            
+            UxUp, UxLw = self.__Remesh(Ux[iSec][0], self.stgPt)
+            UyUp, UyLw = self.__Remesh(Uy[iSec][0], self.stgPt)
+            UzUp, UzLw = self.__Remesh(Uz[iSec][0], self.stgPt)
+
+            """plt.plot(xUp, UxUp)
+            plt.plot(xLw, UxLw)
+            plt.plot(self.viscStruct[iSec][1][:,0], Ux[iSec][1])
+            plt.plot(self.invStruct[0][:,0], U_wg[:,0], 'x')
+            plt.title('Ux = f(x)')
+            plt.show()
+            
+            plt.plot(xUp, UyUp)
+            plt.plot(xLw, UyLw)
+            plt.plot(self.viscStruct[iSec][1][:,0], Uy[iSec][1])
+            plt.plot(self.invStruct[0][:,0], U_wg[:,1], 'x')
+            plt.title('Uy = f(x) smooting{:.2e}'.format(self.config['smoothing']))
+            plt.show()
+            
+            plt.plot(xUp, UzUp)
+            plt.plot(xLw, UzLw)
+            plt.plot(self.viscStruct[iSec][1][:,0], Uz[iSec][1])
+            plt.plot(self.invStruct[0][:,0], U_wg[:,2], 'x')
+            plt.title('Uz = f(x)')
+            plt.show()"""
             
-            UxUp, UxLw = self.__Remesh(Ux[iSec][0], Ux[iSec][1], idx, reg)
-            UyUp, UyLw = self.__Remesh(Uy[iSec][0], Uy[iSec][1], idx, reg)
-            UzUp, UzLw = self.__Remesh(Uz[iSec][0], Uz[iSec][1], idx, reg)
+            """plt.plot(xUp, UyUp, 'x')
+            plt.plot(xLw, UyLw, 'o')
+            #plt.plot(self.viscStruct[iSec][1][:,0], Uy[iSec][1])
+            plt.plot(self.invStruct[0][:,0], U_wg[:,1], 'x')
+            #plt.xlim([0.4, 1.2])
+            #plt.ylim([-0.3, 0.2])
+            plt.title('Uy = f(x) smooting{:.2e}'.format(self.config['smoothing']))
+            plt.show()"""
+
 
-            vSolver.SetVelocities(iSec, 0, UxUp, UzUp, UyUp)
-            vSolver.SetVelocities(iSec, 1, UxLw, UzLw, UyLw)
-            vSolver.SetVelocities(iSec, 2, Ux[iSec][2], Uz[iSec][2], Uy[iSec][2])
+            vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
+            vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
+            vSolver.SetVelocities(iSec, 2, Ux[iSec][1], Uy[iSec][1], Uz[iSec][1])
+
+            MeUp, MeLw = self.__Remesh(Me[iSec][0], self.stgPt)
 
-            MeUp, MeLw = self.__Remesh(Me[iSec][0], Me[iSec][1], idx, reg)
             vSolver.SetMe(iSec, 0, MeUp)
             vSolver.SetMe(iSec, 1, MeLw)
-            vSolver.SetMe(iSec, 2, Me[iSec][2])
+            vSolver.SetMe(iSec, 2, Me[iSec][1])
+
+            RhoUp, RhoLw = self.__Remesh(Rho[iSec][0], self.stgPt)
 
-            RhoUp, RhoLw = self.__Remesh(Rho[iSec][0], Rho[iSec][1], idx, reg)
             vSolver.SetRhoe(iSec, 0, RhoUp)
             vSolver.SetRhoe(iSec, 1, RhoLw)
-            vSolver.SetRhoe(iSec, 2, Rho[iSec][2])
+            vSolver.SetRhoe(iSec, 2, Rho[iSec][1])
+
+            DeltaStarPrevUp, DeltaStarPrevLw = self.__Remesh(self.DeltaStarPrev[iSec][0], self.stgPt)
+            vSolver.SetDeltaStarExt(iSec, 0, DeltaStarPrevUp)
+            vSolver.SetDeltaStarExt(iSec, 1, DeltaStarPrevLw)
+            vSolver.SetDeltaStarExt(iSec, 2, self.DeltaStarPrev[iSec][1])
 
-            for iReg in range(3):
-                vSolver.SetxxExt(iSec, iReg, self.xxPrev[iSec][iReg])
-                vSolver.SetDeltaStarExt(iSec, iReg, self.DeltaStarPrev[iSec][iReg])
-        self.tms['InvToVisc'].stop()
+            xxPrevUp, xxPrevLw = self.__Remesh(self.xxPrev[iSec][0], self.stgPt)
+            vSolver.SetxxExt(iSec, 0, xxPrevUp)
+            vSolver.SetxxExt(iSec, 1, xxPrevLw)
+            vSolver.SetxxExt(iSec, 2, self.xxPrev[iSec][1])
+            
+            
 
     def SetBlowingVel(self, vSolver):
-        self.tms['ViscToInv'].start()
-        BlowingVel = [[np.zeros(0) for iReg in range(len(self.Sections[iSec]))] for iSec in range(len(self.Sections))]
-        xCoord = [[np.zeros(0) for iReg in range(len(self.Sections[iSec]))] for iSec in range(len(self.Sections))]
-        for iSec in range(len(self.Sections)):
-            for iRegion in range(len(self.Sections[iSec])):
-                xCoord[iSec][iRegion] = np.asarray(vSolver.ExtractxCoord(iSec, iRegion))
+        BlowingVel = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
+        for iSec in range(len(self.viscStruct)):
+            for iRegion in range(3):
                 BlowingVel[iSec][iRegion] = np.asarray(vSolver.ExtractBlowingVel(iSec, iRegion))
-                self.DeltaStarPrev[iSec][iRegion] = vSolver.ExtractDeltaStar(iSec, iRegion)
-                self.xxPrev[iSec][iRegion] = vSolver.Extractxx(iSec, iRegion)
-    
-        if self.nDim == 2:
-
-            # Compute cell centroid x positions (in the viscous mesh)
-
-            xCellsUp = np.zeros(len(xCoord[0][0])-1)
-            xCellsUp[0] = (xCoord[0][0][1] - xCoord[0][0][0])/2
-            cumDx = xCoord[0][0][1] - xCoord[0][0][0]
-            for iCell in range(1, len(xCoord[0][0])-1):
-                xCellsUp[iCell] = cumDx + (xCoord[0][0][iCell+1] - xCoord[0][0][iCell])/2
-                cumDx += (xCoord[0][0][iCell] - xCoord[0][0][iCell-1]) # Cumulative space between cells
-
-            xCellsLw = np.zeros(len(xCoord[0][1])-1)
-            xCellsLw[0] = (xCoord[0][1][1] - xCoord[0][1][0])/2
-            cumDx = xCoord[0][1][1] - xCoord[0][0][0]
-            for iCell in range(1, len(xCoord[0][1])-1):
-                xCellsLw[iCell] = cumDx + (xCoord[0][1][iCell+1] - xCoord[0][1][iCell])/2
-                cumDx += (xCoord[0][1][iCell] - xCoord[0][1][iCell-1]) # Cumulative space between cells
-
-            xCellsWk = np.zeros(len(xCoord[0][2])-1)
-            xCellsWk[0] = self.Sections[iSec][0][-1, 0] + (xCoord[0][2][1] - xCoord[0][2][0])/2
-            cumDx = self.Sections[iSec][0][-1, 0] + xCoord[0][2][1] - xCoord[0][2][0]
-            for iCell in range(1, len(xCoord[0][2])-1):
-                xCellsWk[iCell] = cumDx + (xCoord[0][2][iCell+1] - xCoord[0][2][iCell])/2
-                cumDx += (xCoord[0][2][iCell] - xCoord[0][2][iCell-1]) # Cumulative space between cells
-
-            # Create interpolation functions
-
-            blwSplnUp = interp1d(xCellsUp, BlowingVel[0][0], bounds_error=False, fill_value="extrapolate")
-            blwSplnLw = interp1d(xCellsLw, BlowingVel[0][1], bounds_error=False, fill_value="extrapolate")
-            blwSplnWk = interp1d(xCellsWk, BlowingVel[0][2], bounds_error=False, fill_value="extrapolate")
-            
-            # Airfoil (Cell centroids are now computed in the inviscid mesh)
-            for iElm, e in enumerate(self.bndWing.tag.elems):
 
-                # Compute cg of element e
-                posCg = np.zeros(3)
+        if self.config['nDim'] == 2:
+            BlwWing = np.concatenate((BlowingVel[0][0][::-1], BlowingVel[0][1]))
+            BlwWake = BlowingVel[0][2]
+            self.DeltaStarPrev[iSec][0] = np.concatenate((np.asarray(vSolver.ExtractDeltaStar(iSec, 0))[::-1], np.asarray(vSolver.ExtractDeltaStar(iSec, 1))[1:]), axis=None)
+            self.DeltaStarPrev[iSec][1] = np.asarray(vSolver.ExtractDeltaStar(iSec, 2))
+            self.xxPrev[iSec][0] = np.concatenate((np.asarray(vSolver.Extractxx(iSec, 0))[::-1], np.asarray(vSolver.Extractxx(iSec, 1))[1:]), axis=None)
+            self.xxPrev[iSec][1] = np.asarray(vSolver.Extractxx(iSec, 2))
+
+            # Compute cell centers
+            viscCgWing = np.zeros((len(self.viscStruct[0][0][:,0]) - 1, 3))
+            for iElm in range(len(self.viscStruct[0][0][:,0]) - 1):
+                viscCgWing[iElm, 0] = self.viscStruct[0][0][iElm, 0] + self.viscStruct[0][0][iElm + 1, 0]
+                viscCgWing[iElm, 1] = self.viscStruct[0][0][iElm, 1] + self.viscStruct[0][0][iElm + 1, 1]
+                viscCgWing[iElm, 2] = self.viscStruct[0][0][iElm, 2] + self.viscStruct[0][0][iElm + 1, 2]
+            viscCgWing/=2
+
+            viscCgWake = np.zeros((len(self.viscStruct[0][1][:,0]) - 1, 3))
+            for iElm in range(len(self.viscStruct[0][1][:,0]) - 1):
+                viscCgWake[iElm, 0] = self.viscStruct[0][1][iElm, 0] + self.viscStruct[0][1][iElm + 1, 0]
+                viscCgWake[iElm, 1] = self.viscStruct[0][1][iElm, 1] + self.viscStruct[0][1][iElm + 1, 1]
+                viscCgWake[iElm, 2] = self.viscStruct[0][1][iElm, 2] + self.viscStruct[0][1][iElm + 1, 2]
+            viscCgWake/=2
+
+            invCgWing = np.zeros((len(self.bndWing.tag.elems), 3))
+            for iElm, e in enumerate(self.bndWing.tag.elems):
                 for n in e.nodes:
-                    posCg[0] += n.pos[0]
-                    posCg[1] += n.pos[1]
-                    posCg[2] += n.pos[2]
-                posCg/=len(e.nodes)
-                
-                if posCg[1]>=0:
-                    self.bndWing.setU(iElm, float(blwSplnUp(posCg[0])))
-                else:
-                    self.bndWing.setU(iElm, float(blwSplnLw(posCg[0])))
-            
+                    invCgWing[iElm, 0] +=n.pos[0]
+                    invCgWing[iElm, 1] +=n.pos[1]
+                    invCgWing[iElm, 2] +=n.pos[2]
+                invCgWing[iElm, :] /= e.nodes.size()
+            invCgWake = np.zeros((len(self.bndWake.tag.elems), 3))
             for iElm, e in enumerate(self.bndWake.tag.elems):
-
-                # Compute cg of element e
-                posCg = np.zeros(3)
                 for n in e.nodes:
-                    posCg[0] += n.pos[0]
-                    posCg[1] += n.pos[1]
-                    posCg[2] += n.pos[2]
-                posCg/=len(e.nodes)
-
-                self.bndWake.setU(iElm, float(blwSplnWk(posCg[0])))
+                    invCgWake[iElm, 0] +=n.pos[0]
+                    invCgWake[iElm, 1] +=n.pos[1]
+                    invCgWake[iElm, 2] +=n.pos[2]
+                invCgWake[iElm, :] /= e.nodes.size()
+            
+            if self.config['nDim'] == 3:
+                invCgWing[:,[1,2]] = invCgWing[:,[2,1]]
+                invCgWake[:,[1,2]] = invCgWake[:,[2,1]]
             
+            """plt.plot(self.invStruct[0][:,0], self.invStruct[0][:,1], 'x')
+            plt.plot(invCgWing[:, 0], invCgWing[:,1], 'o')
+            plt.show()"""
 
-        elif self.nDim == 3:
-            pass
-        self.tms['ViscToInv'].stop()
+        invCgWingUp = invCgWing[invCgWing[:,1]>0, :]
+        invCgWingLw = invCgWing[invCgWing[:,1]<0, :]
+
+        mappedBlwWing = np.zeros(len(invCgWingUp) + len(invCgWingLw))
+        self.tms['Interpolation'].start()
+        #mappedBlwWing[:self.config['nPoints'][0]] = RBFInterpolator(viscCgWing[:self.config['nPoints'][0],:], BlwWing[:self.config['nPoints'][0]], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=1e-6, degree=self.config['degree'])(invCgWingUp)
+        #mappedBlwWing[self.config['nPoints'][0]:] = RBFInterpolator(viscCgWing[self.config['nPoints'][0]:, :], BlwWing[self.config['nPoints'][0]:], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=1e-6, degree=self.config['degree'])(invCgWingLw)
+        mappedBlwWing[:len(invCgWingUp)] = RBFInterpolator(viscCgWing[:self.config['nPoints'][0],:], BlwWing[:self.config['nPoints'][0]], neighbors=1, kernel='linear', smoothing=1e-6, degree=0)(invCgWingUp)
+        mappedBlwWing[len(invCgWingUp):] = RBFInterpolator(viscCgWing[self.config['nPoints'][0]:, :], BlwWing[self.config['nPoints'][0]:], neighbors=1, kernel='linear', smoothing=1e-6, degree=0)(invCgWingLw)
+        #mappedBlwWake = RBFInterpolator(viscCgWake, BlwWake, neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(invCgWake)
+        mappedBlwWake = RBFInterpolator(viscCgWake, BlwWake, neighbors=1, kernel='linear', smoothing=1e-6, degree=0)(invCgWake)
+
+        """plt.plot(invCgWing[:,0], invCgWing[:,1])
+        plt.show()"""
+        """plt.plot(invCgWing[:, 0], mappedBlwWing, 'o')
+        plt.plot(viscCgWing[:,0], BlwWing, 'x')
+        #plt.xlim([0.8, 1.2])
+        plt.show()"""
+
+        """interpWing = NearestNDInterpolator(list(zip(viscCgWing[:,0], viscCgWing[:,1], viscCgWing[:,2])), BlwWing)
+        mappedBlwWing = interpWing(invCgWing[:,0], invCgWing[:,1], invCgWing[:,2])
+        interpWake = NearestNDInterpolator(list(zip(viscCgWake[:,0], viscCgWake[:,1], viscCgWake[:,2])), BlwWake)
+        mappedBlwWake = interpWake(invCgWake[:,0], invCgWake[:,1], invCgWake[:,2])"""
+
+        self.tms['Interpolation'].stop()
+        for iElm in range(len(self.bndWing.tag.elems)):
+            self.bndWing.setU(iElm, mappedBlwWing[iElm])
+        for iElm in range(len(self.bndWake.tag.elems)):
+            self.bndWake.setU(iElm, mappedBlwWake[iElm])
 
 
     def RunSolver(self):
@@ -224,227 +264,266 @@ class Interpolator:
     def GetCm(self):
         return self.iSolver.Cm
 
-    def __RealToSections(self, var, varWk):
-        pass
+    def __Remesh(self, vec, idx_stag):
+        xUp = vec[:idx_stag+1]
+        xUp = xUp[::-1]
 
-    def __InterpToSections(self, var, varWk):
+        xLw = vec[idx_stag:]
+        return xUp, xLw
 
+    def __RbfToSections(self, var, varWk):
         varWk = np.reshape(varWk, [len(varWk), 1])
-        completeWake = np.hstack((self.Wake, varWk))
+        completeWake = np.hstack((self.invStruct[1], varWk))
         _, idx = np.unique(completeWake[:,:3], axis=0, return_index=True)
+
+        saveWake = np.where(completeWake[:,0]==1)
+        saveWake = saveWake[0]
         completeWake = completeWake[np.sort(idx)]
 
         var = np.reshape(var, [len(var), 1])
-        completeWing = np.hstack((self.Wing, var))
+        completeWing = np.hstack((self.invStruct[0], var))
 
-        wingUp = completeWing[completeWing[:,2]>=0, :]
-        wingLw = completeWing[completeWing[:,2]<=0, :]
+        delIdx = np.where(completeWing[:,0]==1.0)
+        delIdx = delIdx[0]
+        save = completeWing[delIdx, 3]
+        print('save', save)
 
-        mappedVar = [[np.zeros(0) for _ in range(3)] for _ in range(len(self.Sections))]
+        completeWing = np.delete(completeWing, delIdx[0], axis = 0)
+        completeWingUp = completeWing[completeWing[:,1]>0,:]
+        completeWingLw = completeWing[completeWing[:,1]<0,:]
 
-        for iSec in range(len(mappedVar)):
-            for iReg in range(len(mappedVar[iSec])):
-                if iReg==0:
-                    interpFunc = interp1d(wingUp[:,0], wingUp[:,3])
-                if iReg==1:
-                    interpFunc = interp1d(wingLw[:,0], wingLw[:,3])
-                if iReg==2:
-                    interpFunc = interp1d(completeWake[:,0], completeWake[:,3])
-                
-                mappedVar[iSec][iReg] = interpFunc(self.Sections[iSec][iReg][:,0])
-
-        return mappedVar
-
-
-    def __RbfToSections(self, var, varWk):
-
-        varWk = np.reshape(varWk, [len(varWk), 1])
-        completeWake = np.hstack((self.Wake, varWk))
-        _, idx = np.unique(completeWake[:,:3], axis=0, return_index=True)
-        completeWake = completeWake[np.sort(idx)]
+        completeWingUp = np.vstack((completeWingUp, completeWing[delIdx[0], :]))
+        completeWingLw = np.vstack((completeWingLw, completeWing[delIdx[0], :]))
 
-        var = np.reshape(var, [len(var), 1])
-        completeWing = np.hstack((self.Wing, var))
-
-        wingUp = completeWing[completeWing[:,2]>=0]
-        wingLw = completeWing[completeWing[:,2]<=0]
-
-        mappedVar = [[np.zeros(0) for _ in range(3)] for _ in range(len(self.Sections))]
+        mappedVar = [[np.zeros(0) for _ in range(3)] for _ in range(len(self.viscStruct))]
 
         for iSec in range(len(mappedVar)):
             for iReg in range(len(mappedVar[iSec])):
-                if iReg==0:
-                    mappedVar[iSec][iReg] = RBFInterpolator(wingUp[:,:3], wingUp[:,3], neighbors=10, smoothing=0.1, epsilon=0.6, kernel='gaussian', degree=1)(self.Sections[iSec][iReg])
-                    #mappedVar[iSec][iReg] = RBFInterpolator(wingUp[:,:3], wingUp[:,3])(self.Sections[iSec][iReg])
-                if iReg==1:
-                    mappedVar[iSec][iReg] = RBFInterpolator(wingLw[:,:3], wingLw[:,3], neighbors=10, smoothing=0.1, epsilon=0.6, kernel='gaussian', degree=1)(self.Sections[iSec][iReg])
-                    #mappedVar[iSec][iReg] = RBFInterpolator(wingLw[:,:3], wingLw[:,3])(self.Sections[iSec][iReg])
-                if iReg==2:
-                    mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:2], completeWake[:,3], neighbors=10, smoothing=0.1, kernel='linear', degree=0)(self.Sections[iSec][iReg][:,:2])
-        return mappedVar
-
-
-    def __Remesh(self, vectorUp, vectorLw, idx, reg):
-        """Rebuilds input vector to match the stagnation point position.
-        
-        args:
-        ----
-        - VectorUp/VectorLw upper/lower side vector of the desired variable from x_glob = 0 to x_glob = 1
-        - idx : Stagnation point position from leading edge point (x_glob=0)
-        - reg : 0/1 indicates on which side the stagnation point is.
-
-        return : vUp/vLw distribution from x_glob=stag to x_glob=1 on upper/lower sides.
-        ------
-        """
-
-        if reg==0:
-            if idx!=0:
-                vLw = np.insert(vectorLw, 0, np.flip(vectorUp[1:idx]))
-                vUp = np.delete(vectorUp, np.s_[:idx-1])
-            else:
-                vLw = vectorLw
-                vUp = vectorUp
-        elif reg==1:
-            if idx!=0:
-                vUp = np.insert(vectorUp, 0, np.flip(vectorLw[1:idx]))
-                vLw = np.delete(vectorLw, np.s_[:idx-1])
-            else:
-                vUp = vectorUp
-                vLw = vectorLw
-        else:
-            raise RuntimeError('Incorrect region specified during remeshing.')
-        
-        return vUp, vLw
+                if iReg==0: # Wing
+                    mappedVar[iSec][iReg] = np.zeros(len(self.viscStruct[iSec][0][:,0]))
+                    mappedVar[iSec][iReg][:self.config['nPoints'][0]] = RBFInterpolator(completeWingUp[:,:2], completeWingUp[:,3], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(self.viscStruct[iSec][0][:self.config['nPoints'][0],:2])
+                    mappedVar[iSec][iReg][self.config['nPoints'][0]:] = RBFInterpolator(completeWingLw[:,:2], completeWingLw[:,3], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(self.viscStruct[iSec][0][self.config['nPoints'][0]:,:2])
+                    
+                    #interpWWW = NearestNDInterpolator(list(zip(completeWing[:,0], completeWing[:,1], completeWing[:,2])), completeWing[:,3])
+                    #mappedVar[iSec][iReg] = interpWWW(self.viscStruct[iSec][0][:,0], self.viscStruct[iSec][0][:,1], self.viscStruct[iSec][0][:,2])
+                    
+                    #interp = CloughTocher2DInterpolator(list(zip(completeWing[:,0], completeWing[:,1])), completeWing[:,3])
+                    #mappedVar[iSec][iReg] = interp(self.viscStruct[0][0][:,0], self.viscStruct[0][0][:,1])
+
+                    #mappedVar[iSec][iReg] = griddata(completeWing[:,:2], completeWing[:,3], self.viscStruct[iSec][0][:,:2])
+
+                if iReg==1: # Wake
+                    
+                    #mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:3], completeWake[:,3], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(self.viscStruct[iSec][1])
+                    mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:3], completeWake[:,3], neighbors=1, kernel='linear', smoothing=1e-6, degree=0)(self.viscStruct[iSec][1])
+                    
+                    #interpWWW = NearestNDInterpolator(list(zip(completeWake[:,0], completeWake[:,1], completeWake[:,2])), completeWake[:,3])
+                    #mappedVar[iSec][iReg] = interpWWW(self.viscStruct[iSec][1][:,0], self.viscStruct[iSec][1][:,1], self.viscStruct[iSec][1][:,2])
         
+                    #mappedVar[iSec][iReg] = griddata(completeWake[:,:2], completeWake[:,3], self.viscStruct[iSec][1][:,:2], method='linear')
+                    pass
 
-    def __GetWing(self, nSections, nDim, nPoints, eps, k):
-        "Extract and sort the node coordinates of the airfoil."
-
-        if nDim == 2 and nSections != 1:
-            raise RuntimeError('Cannot create more than 1 section on 2D geometry. Specified:', nSections)
-            
-        wing = np.zeros((len(self.bndWing.nodes),3))
-        wake = np.zeros((len(self.bndWake.nodes),3))
 
-        # Extract coordinates.
-        for iNode, n in enumerate(self.bndWing.nodes):
-            for iDir in range(len(wing[0,:])):
-                wing[iNode,iDir] = n.pos[iDir]
+        #mappedVar[iSec][0][0] = save[0]
+        #mappedVar[iSec][0][-1] = save[-1]
+        #mappedVar[iSec][1][0] = saveWake[0]
+        """fig = plt.figure()
+        ax = fig.add_subplot(projection='3d')
+        ax.scatter(completeWing[:,0], completeWing[:,1], completeWing[:,2], completeWing[:,3], s=0.2)
+        ax.set_ylim([-0.5, 0.5])
+        ax.set_xlabel('x')
+        ax.set_ylabel('y')
+        ax.set_zlabel('z')
+        plt.show()"""
+        #plt.plot(completeWing[:,0], completeWing[:,3], 'x')
+        #plt.plot(completeWake[:,0], completeWake[:,3], 'x')
+        """plt.plot(self.viscStruct[0][0][:,0], mappedVar[0][0])
+        plt.plot(self.viscStruct[0][1][:,0], mappedVar[0][1])
+        plt.show()"""
+        return mappedVar
 
-        for iNode, n in enumerate(self.bndWake.nodes):
-            for iDir in range(len(wake[0,:])):
-                wake[iNode, iDir] = n.pos[iDir]
 
-        if nDim == 2: # 2D case
-            wing[:,[1,2]] = wing[:,[2,1]]
-            wake[:,[1,2]] = wake[:,[2,1]]
-            
-        # Separate upper and lower sides
-
-        wingUp = wing[wing[:,2]>=0]
-        wingLw = wing[wing[:,2]<=0]
-
-        # Sort upper/lower side and remove duplicates
-        wingUp = np.delete(wingUp, wingUp[:,1]>1, 0)
-        wingLw = np.delete(wingLw, wingLw[:,1]>1, 0)
-
-        _, idxUp = np.unique(wingUp, axis=0, return_index=True)
-        wingUp = wingUp[np.sort(idxUp)]
-        _, idxLw = np.unique(wingLw, axis=0, return_index=True)
-        wingLw = wingLw[np.sort(idxLw)]
-
-        if nDim == 2:
-            arfSplnUp = interp1d(wingUp[:,0], wingUp[:,2], kind='cubic')
-            arfSplnLw = interp1d(wingLw[:,0], wingLw[:,2], kind='cubic')
-
-            Sections = [[np.zeros((nPoints[1] if i==2 else nPoints[0], 3)) for i in range(3)]]
-            for iReg in range(len(Sections[0])):
-                Sections[0][iReg][:,1] = wingUp[0,1]
-                if iReg != 2:
-                    Sections[0][iReg][:,0] = np.linspace(min(wing[:,0]), max(wing[:,0]), nPoints[0])
-
-                    if iReg == 0:
-                        for iPoint, x in enumerate(Sections[0][iReg][:,0]):
-                            Sections[0][iReg][iPoint,2] = arfSplnUp(x)
-                    elif iReg == 1:
-                        for iPoint, x in enumerate(Sections[0][iReg][:,0]):
-                            Sections[0][iReg][iPoint,2] = arfSplnLw(x)
-                else:
-                    Sections[0][iReg][:,0] = np.linspace(min(wake[:,0]), max(wake[:,0]), nPoints[1])
+    def __GetWing(self, bnd_wing, bnd_wake, config):
 
-        
-        elif nDim == 3:
-            sUp = (len(wingUp[:,0])-np.sqrt(2*len(wingUp[:,0])) + len(wingUp[:,0])+np.sqrt(2*len(wingUp[:,0])))/2
-            sLw = (len(wingLw[:,0])-np.sqrt(2*len(wingLw[:,0])) + len(wingLw[:,0])+np.sqrt(2*len(wingLw[:,0])))/2
+        if config['nDim'] == 2 and config['nSections'] != 1:
+            raise RuntimeError('Cannot create more than 1 section on 2D geometry. Specified:', config['nSections'])
 
-            tckUp = bisplrep(wingUp[:,0], wingUp[:,1], wingUp[:,2], kx=k[0], ky=k[1], s=sUp)
-            tckLw = bisplrep(wingLw[:,0], wingLw[:,1], wingLw[:,2], kx=k[0], ky=k[1], s=sLw)
+        viscStruct = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0], 3)) for i in range(2)] for _ in range(config['nSections'])]
 
-            zDistri = np.linspace(min(wing[:,1]), max(wing[:,1]), nSections)
+        invStruct = [np.zeros((len(bnd_wake.nodes) if iReg==1 else len(bnd_wing.nodes), 3)) for iReg in range(2)]
 
-            zDistri[-1] -= 0.2
-            
-            Sections = [[np.zeros((nPoints[1] if i==2 else nPoints[0], 3)) for i in range(3)] for _ in range(nSections)]
-
-            for iSec, y in enumerate(zDistri):
-                for iReg in range(3):
-                    print('Sec', iSec, 'y', y, min(wing[abs(wing[:,1]-y)<eps, 0]), max(wing[abs(wing[:,1]-y)<eps, 0]))
-                    print(wing[abs(wing[:,1]-y)<0.02, :])
-                    if iReg == 0:
-                        Sections[iSec][iReg][:,0] = np.linspace(min(wing[abs(wing[:,1]-y)<eps, 0]), max(wing[abs(wing[:,1]-y)<eps, 0]), nPoints[0])
-                    elif iReg == 1:
-                        Sections[iSec][iReg][:,0] = np.linspace(min(wing[abs(wing[:,1]-y)<eps, 0]), max(wing[abs(wing[:,1]-y)<eps, 0]), nPoints[0])
-                    elif iReg == 2:
-                        #print('sec', iSec, 'min', min(wake[:, 0]), 'max', max(wake[:, 0]))
-                        Sections[iSec][iReg][:,0] = np.linspace(max(wing[abs(wing[:,1]-y)<eps, 0]), max(wake[:, 0]), nPoints[1])
-                    Sections[iSec][iReg][:,1] = y
-                    for iPt in range(len(Sections[iSec][iReg][:,0])):
-                        if iReg==0:
-                            Sections[iSec][iReg][iPt,2] = bisplev(Sections[iSec][iReg][iPt, 0], Sections[iSec][iReg][iPt, 1], tckUp)
-                        elif iReg==1:
-                            Sections[iSec][iReg][iPt,2] = bisplev(Sections[iSec][iReg][iPt, 0], Sections[iSec][iReg][iPt, 1], tckLw)
-                        elif iReg == 2:
-                            Sections[iSec][iReg][iPt,2] = 0
-                        else:
-                            raise RuntimeError("dart::vInterpolator3 l.83 0<=iReg<=2")
-
-        # Add leading edge point
-        #Sections.pop()
-
-        for iSec in range(len(Sections)):
-            for iReg in range(2):
-                Sections[iSec][iReg][0,2] = 0
+        for iNode, n in enumerate(bnd_wing.nodes):
+            invStruct[0][iNode, 0] = n.pos[0]
+            invStruct[0][iNode, 1] = n.pos[1]
+            invStruct[0][iNode, 2] = n.pos[2]
+        for iNode, n in enumerate(bnd_wake.nodes):
+            invStruct[1][iNode, 0] = n.pos[0]
+            invStruct[1][iNode, 1] = n.pos[2]
+            invStruct[1][iNode, 2] = n.pos[1]
         
-        """fig = plt.figure()
+        """xy = invStruct[0][abs(invStruct[0][:,1]-0.5) <= 5e-2, :]
+        fig = plt.figure()
         ax = fig.add_subplot(projection='3d')
-        ax.scatter(wing[:,0], wing[:,1], wing[:,2], s=0.3)
-        ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)
-        for iSec in range(len(Sections)):
-            ax.scatter(Sections[iSec][0][:,0], Sections[iSec][0][:,1], Sections[iSec][0][:,2], color='red')
-            ax.scatter(Sections[iSec][1][:,0], Sections[iSec][1][:,1], Sections[iSec][1][:,2], color='red')
-            ax.scatter(Sections[iSec][2][:,0], Sections[iSec][2][:,1], Sections[iSec][2][:,2], color='red')
-        ax.set_zlim([-0.5, 0.5])
+        ax.scatter(xy[:,0], xy[:,1], xy[:,2], s=0.2)
+        ax.set_ylim([-0.5, 0.5])
         ax.set_xlabel('x')
         ax.set_ylabel('y')
         ax.set_zlabel('z')
         plt.show()
 
-        plt.plot(Sections[0][0][:,0], Sections[0][0][:,2], 'x-', color='red')
-        plt.plot(Sections[0][1][:,0], Sections[0][1][:,2], 'x-', color='red')
-        plt.plot(wingUp[:,0], wingUp[:,2], 'o', color='blue')
-        plt.plot(wingLw[:,0], wingLw[:,2], 'o', color='blue')
-        plt.ylim([-0.5, 0.5])
-        plt.show()
-        #quit()
+        plt.plot(xy[:,0], xy[:,2], 'x')
+        plt.show()"""
+
+        if config['nDim'] == 2:
+            spanDistri = np.zeros(1)
+            spanDistri[0] = invStruct[0][0,2] # Span direction is the z direction 
+
+        elif config['nDim'] == 3:
+            invStruct[0][:,[1,2]] = invStruct[0][:,[2,1]] # Inverse y and z so that spanwise direction is in z
+            invStruct[1][:,[1,2]] = invStruct[1][:,[2,1]] # Inverse y and z so that spanwise direction is in z
+            spanDistri = np.linspace(min(invStruct[0][:,2]), max(invStruct[0][:,2]), config['nSections'])
+            spanDistri[0] +=0.1
+            spanDistri[-1] -=0.1
+
+        for iSec, z in enumerate(spanDistri):
+
+            # Create x distribution at the corresponding span position
+            minX = min(invStruct[0][abs(invStruct[0][:,2]-z)<=config['eps'], 0])
+            chordLength = max(invStruct[0][abs(invStruct[0][:,2]-z)<=config['eps'], 0]) - minX
+            xDistri = minX + 1/2*(1 + np.cos(np.linspace(0, np.pi, config['nPoints'][0])))*chordLength
+            nElmSide = len(xDistri)
+            viscStruct[iSec][0][:,0] = np.concatenate((xDistri, xDistri[::-1]), axis=None)
+            viscStruct[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
+
+            wakeLength = np.max(invStruct[1][:,0]) - np.max(viscStruct[iSec][0][:,0])
+            viscStruct[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
+            viscStruct[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
+
+            if config['nDim'] == 2:
+                upper, lower = self.__SeparateArf(invStruct[0][:,0], invStruct[0][:,1])
+                viscStruct[iSec][0][:nElmSide,1] = interp1d(upper[:,0], upper[:,1])(viscStruct[iSec][0][:nElmSide,0])
+                viscStruct[iSec][0][nElmSide:,1] = interp1d(lower[:,0], lower[:,1])(viscStruct[iSec][0][nElmSide:,0])
+                viscStruct[iSec][0] = np.delete(viscStruct[iSec][0], nElmSide, axis=0)
+
+                viscStruct[iSec][1][:,1] = interp1d(invStruct[1][:,0], invStruct[1][:,1])(viscStruct[iSec][1][:,0])
+
+                """plt.plot(viscStruct[iSec][0][:,0], viscStruct[iSec][0][:,1], '-')
+                plt.plot(viscStruct[iSec][0][:nElmSide,0], viscStruct[iSec][0][:nElmSide,1], 'x')
+                plt.plot(viscStruct[iSec][0][nElmSide:,0], viscStruct[iSec][0][nElmSide:,1], 'o')
+                plt.plot(viscStruct[iSec][1][:,0], np.zeros(len(viscStruct[iSec][1][:,0])), 'x-')
+                plt.xlim([-0.5,2])
+                plt.ylim([-0.5, 0.5])
+                plt.show()"""
+
+            elif config['nDim'] == 3:
+                nPointsMin = (config['k'][0]+1) * (config['k'][1]+1) # See Scipy bisplrep doc
+
+                wingSpan = max(invStruct[0][:,2]) - min(invStruct[0][:,2])
+                portion = invStruct[0][abs(invStruct[0][:,2] - z) < 0.01 * wingSpan,:]
+
+                portionUpper, portionLower = self.__SeparateArf(portion[:,0], portion[:,1], portion[:,2])
+                if len(portionUpper)<=nPointsMin or len(portionLower)<=nPointsMin:
+                    print(ccolors.ANSI_YELLOW,'dart::Interpolator Warning: Found {:.0f} points within 1% of spanwise position z={:.2f}. Trying 10%. '.format(
+                        len(portion), z),ccolors.ANSI_RESET)
+                    portion = invStruct[0][abs(invStruct[0][:,2] - z) < 0.1 * wingSpan,:]
+                    if len(portion)<=nPointsMin:
+                        print(ccolors.ANSI_RED, 'Failed. Continuing\n', ccolors.ANSI_RESET)
+                    else:
+                        print(ccolors.ANSI_GREEN, 'Sucess.\n', ccolors.ANSI_RESET)
+
+                # Upper
+                sWing = (len(portionUpper[:,0])-np.sqrt(2*len(portionUpper[:,0])) + len(portionUpper[:,0])+np.sqrt(2*len(portionUpper[:,0])))/2
+                tckWingUp = bisplrep(portionUpper[:,0], portionUpper[:,2], portionUpper[:,1], kx=config['k'][0], ky=config['k'][1], s=sWing)
+                for iPoint in range(nElmSide):
+                    viscStruct[iSec][0][iPoint,1] = bisplev(viscStruct[iSec][0][iPoint,0], viscStruct[iSec][0][iPoint,2], tckWingUp)
+                
+                # Lower
+                sWing = (len(portionLower[:,0])-np.sqrt(2*len(portionLower[:,0])) + len(portionLower[:,0])+np.sqrt(2*len(portionLower[:,0])))/2
+                tckWingUp = bisplrep(portionLower[:,0], portionLower[:,2], portionLower[:,1], kx=config['k'][0], ky=config['k'][1], s=sWing)
+                for iPoint in range(nElmSide, len(viscStruct[iSec][0][:,0])):
+                    viscStruct[iSec][0][iPoint,1] = bisplev(viscStruct[iSec][0][iPoint,0], viscStruct[iSec][0][iPoint,2], tckWingUp)
+                
+                """invStruct[1] = invStruct[1][invStruct[1][:, 0].argsort()]
+                fig = plt.figure()
+                ax = fig.add_subplot(projection='3d')
+                ax.scatter(invStruct[1][:,0], invStruct[1][:,1], invStruct[1][:,2], s=0.2)
+                ax.set_ylim([-0.5, 0.5])
+                ax.set_xlabel('x')
+                ax.set_ylabel('y')
+                ax.set_zlabel('z')
+                plt.show()"""
+                # Wake
+                """sWake = (len(invStruct[1][:,0])-np.sqrt(2*len(invStruct[1][:,0])) + len(invStruct[1][:,0])+np.sqrt(2*len(invStruct[1][:,0])))/2
+                tckWake = bisplrep(invStruct[1][:,0], invStruct[1][:,2], invStruct[1][:,1], kx=config['k'][0], ky=config['k'][1], s=sWake)"""
+                for iPoint in range(len(viscStruct[iSec][1][:,0])):
+                    #viscStruct[iSec][1][iPoint,1] = bisplev(viscStruct[iSec][1][iPoint,0], viscStruct[iSec][1][iPoint,2], tckWake)
+                    viscStruct[iSec][1][iPoint,1] = 0
+
+        if config['nDim'] ==3:
+            fig = plt.figure()
+            ax = fig.add_subplot(projection='3d')
+            ax.scatter(invStruct[0][:,0], invStruct[0][:,1], invStruct[0][:,2], s=0.2)
+            for iSec in range(len(viscStruct)):
+                ax.scatter(viscStruct[iSec][0][:,0], viscStruct[iSec][0][:,1], viscStruct[iSec][0][:,2], s=0.9, color='red')
+            ax.set_ylim([-0.5, 0.5])
+            ax.set_xlabel('x')
+            ax.set_ylabel('y')
+            ax.set_zlabel('z')
+            plt.show()
+        return invStruct, viscStruct
+
+    def __SeparateArf(self, xInp, yInp, zInp=None):
+        if zInp is None:
+            zInp = np.zeros(len(xInp))
+        if len(xInp) != len(yInp) or len(xInp) != len(zInp):
+            raise RuntimeError ('dart::Interpolator::__SeparateArf Missmatch between intput vector sizes.')
+
+        Q=np.zeros((len(xInp), 3))
+        Q[:,0] = xInp.copy()
+        Q[:,1] = yInp.copy()
+        Q[:,2] = zInp.copy()
+
+        Q = Q[Q[:, 0].argsort()]
+
+        le = Q[0,:]
+        te = Q[-1,:]
+        ymean = (le[1] + te[1]) / 2
+
+        upper, lower, save = le, le, le
+
+        flag1 = 0
+        iPoint = 0
+        for i in range(len(Q[:,0])):
+            if Q[iPoint,1] == 0.0:
+                if flag1 == 1:
+                    Q=np.delete(Q, iPoint, axis = 0)
+                    iPoint-=1
+                else:
+                    flag1 = 1
+            iPoint +=1
+
+        for iPoint in range(1, len(Q[:,0])):
+            if Q[iPoint,1] > ymean:
+                upper = np.vstack((upper, Q[iPoint,:]))
+            elif Q[iPoint, 1] < ymean:
+                lower = np.vstack((lower, Q[iPoint,:]))
+            else:
+                save = np.vstack((save, Q[iPoint, :]))
+        
+        _, idx = np.unique(save, return_index=True, axis=0)
+        save=save[np.sort(idx)]
 
-        for iSec in range(len(Sections)):
-            plt.plot(Sections[iSec][0][:,0], Sections[iSec][0][:,2], 'x-', color='red')
-            plt.plot(Sections[iSec][1][:,0], Sections[iSec][1][:,2], 'x-', color='red')
-            plt.xlim([-0.5, 1.5])
-            plt.ylim([-0.5,0.5])
-            plt.show()"""
-        return Sections, wing, wake
+        #save = np.delete(save, 0, axis=0)
 
+        
+        try:
+            if len(save[:,0])>0:
+                upper = np.vstack((upper, save))
+                lower = np.vstack((lower, save))
+        except:
+            pass
 
+        upper = np.vstack((upper, te))
+        lower = np.vstack((lower, te))
 
-    
\ No newline at end of file
+        return upper, lower
diff --git a/dart/pyVII/vInterpolator2.py b/dart/pyVII/vInterpolator2.py
deleted file mode 100644
index ab546c8ee6883727cc0e4281d37eda712b643e9f..0000000000000000000000000000000000000000
--- a/dart/pyVII/vInterpolator2.py
+++ /dev/null
@@ -1,396 +0,0 @@
-import numpy as np
-from matplotlib import pyplot as plt
-from fwk.coloring import ccolors
-
-# 3D spline with Bézier curves
-from scipy.interpolate import bisplrep
-from scipy.interpolate import bisplev
-
-# 2D spline with Bézier curves
-from scipy.interpolate import splprep
-from scipy.interpolate import splev
-from scipy.integrate import cumtrapz
-
-from scipy.interpolate import interp1d
-
-from scipy.interpolate import RBFInterpolator
-
-import fwk
-
-
-
-class Interpolator:
-    def __init__(self, _dartSolver, bnd_wing, bnd_wake, _cfg):
-        self.iSolver = _dartSolver
-
-        self.bndWing = bnd_wing
-        self.bndWake = bnd_wake
-        self.invStruct, self.viscStruct = self.__GetWing(bnd_wing, bnd_wake, _cfg)
-        self.config = _cfg
-
-        self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
-        self.DeltaStarPrev = [[np.zeros(len(self.viscStruct[iSec][iReg][:,0])) for iReg in range(2)] for iSec in range(len(self.viscStruct))]
-        self.tms = fwk.Timers()
-
-    def GetInvBC(self, vSolver, couplIter):
-        self.tms['InvToVisc'].start()
-
-        # Extract parameters from the inviscid solver
-        # Particular care must be taken with space directions; In the inviscid, the airfoil is in the x,z plane.
-        # While setting up the viscous solver, y and z must be inverted. 
-
-        # Wing
-
-        nElm = len(self.bndWing.nodes)
-        U_wg, M_wg, Rho_wg = np.zeros((nElm, 3)), np.zeros(nElm), np.zeros(nElm)
-        for iNode, n in enumerate(self.bndWing.nodes):
-            U_wg[iNode,0] = self.iSolver.U[n.row][0]
-            U_wg[iNode,1] = self.iSolver.U[n.row][1]
-            U_wg[iNode,2] = self.iSolver.U[n.row][2]
-            M_wg[iNode] = self.iSolver.M[n.row]
-            Rho_wg[iNode] = self.iSolver.rho[n.row]
-        
-        # Wake
-
-        nElm = len(self.bndWake.nodes)
-        U_wk, M_wk, Rho_wk = np.zeros((nElm, 3)), np.zeros(nElm), np.zeros(nElm)
-        for iNode, n in enumerate(self.bndWake.nodes):
-            U_wk[iNode,0] = self.iSolver.U[n.row][0]
-            U_wk[iNode,1] = self.iSolver.U[n.row][2] # x and y inverted
-            U_wk[iNode,2] = self.iSolver.U[n.row][1]
-            M_wk[iNode] = self.iSolver.M[n.row]
-            Rho_wk[iNode] = self.iSolver.rho[n.row]
-
-        Ux = self.__RbfToSections(U_wg[:,0], U_wk[:,0])
-        Uy = self.__RbfToSections(U_wg[:,1], U_wk[:,1])
-        Uz = self.__RbfToSections(U_wg[:,2], U_wk[:,2])
-
-        Me = self.__RbfToSections(M_wg, M_wk)
-        Rho = self.__RbfToSections(Rho_wg, Rho_wk)
-
-        # Find stagnation point
-        for iSec in range(len(self.viscStruct)):
-
-            #if couplIter == 0:
-            self.stgPt = np.argmin(Ux[iSec][0]**2+Uy[iSec][0]**2)
-
-            xUp, xLw = self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt)
-            yUp, yLw = self.__Remesh(self.viscStruct[iSec][0][:,1], self.stgPt)
-            zUp, zLw = self.__Remesh(self.viscStruct[iSec][0][:,2], self.stgPt)
-
-            vSolver.SetGlobMesh(iSec, 0, xUp, yUp, zUp)
-            vSolver.SetGlobMesh(iSec, 1, xLw, yLw, zLw)
-            vSolver.SetGlobMesh(iSec, 2, self.viscStruct[iSec][1][:,0],
-                                         self.viscStruct[iSec][1][:,1],
-                                         self.viscStruct[iSec][1][:,2])
-            
-            UxUp, UxLw = self.__Remesh(Ux[iSec][0], self.stgPt)
-            UyUp, UyLw = self.__Remesh(Uy[iSec][0], self.stgPt)
-            UzUp, UzLw = self.__Remesh(Uz[iSec][0], self.stgPt)
-
-            vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
-            vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
-            vSolver.SetVelocities(iSec, 2, Ux[iSec][1], Uz[iSec][1], Uy[iSec][1])
-
-            MeUp, MeLw = self.__Remesh(Me[iSec][0], self.stgPt)
-
-            vSolver.SetMe(iSec, 0, MeUp)
-            vSolver.SetMe(iSec, 1, MeLw)
-            vSolver.SetMe(iSec, 2, Me[iSec][1])
-
-            RhoUp, RhoLw = self.__Remesh(Rho[iSec][0], self.stgPt)
-
-            vSolver.SetRhoe(iSec, 0, RhoUp)
-            vSolver.SetRhoe(iSec, 1, RhoLw)
-            vSolver.SetRhoe(iSec, 2, Rho[iSec][1])
-
-            DeltaStarPrevUp, DeltaStarPrevLw = self.__Remesh(self.DeltaStarPrev[iSec][0], self.stgPt)
-            vSolver.SetDeltaStarExt(iSec, 0, DeltaStarPrevUp)
-            vSolver.SetDeltaStarExt(iSec, 1, DeltaStarPrevLw)
-            vSolver.SetDeltaStarExt(iSec, 2, self.DeltaStarPrev[iSec][1])
-
-            xxPrevUp, xxPrevLw = self.__Remesh(self.xxPrev[iSec][0], self.stgPt)
-            vSolver.SetxxExt(iSec, 0, xxPrevUp)
-            vSolver.SetxxExt(iSec, 1, xxPrevLw)
-            vSolver.SetxxExt(iSec, 2, self.xxPrev[iSec][1])
-            
-            
-
-    def SetBlowingVel(self, vSolver):
-        BlowingVel = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
-        for iSec in range(len(self.viscStruct)):
-            for iRegion in range(3):
-                BlowingVel[iSec][iRegion] = np.asarray(vSolver.ExtractBlowingVel(iSec, iRegion))
-
-
-        if self.config['nDim'] == 2:
-            BlwWing = np.concatenate((BlowingVel[0][0][::-1], BlowingVel[0][1]))
-            BlwWake = BlowingVel[0][2]
-            self.DeltaStarPrev[iSec][0] = np.concatenate((np.asarray(vSolver.ExtractDeltaStar(iSec, 0))[::-1], np.asarray(vSolver.ExtractDeltaStar(iSec, 1))[1:]), axis=None)
-            self.DeltaStarPrev[iSec][1] = np.asarray(vSolver.ExtractDeltaStar(iSec, 2))
-            self.xxPrev[iSec][0] = np.concatenate((np.asarray(vSolver.Extractxx(iSec, 0))[::-1], np.asarray(vSolver.Extractxx(iSec, 1))[1:]), axis=None)
-            self.xxPrev[iSec][1] = np.asarray(vSolver.Extractxx(iSec, 2))
-
-            # Compute cell centers
-            viscCgWing = np.zeros((len(self.viscStruct[0][0][:,0]) - 1, 3))
-            for iElm in range(len(self.viscStruct[0][0][:,0]) - 1):
-                viscCgWing[iElm, 0] = self.viscStruct[0][0][iElm, 0] + self.viscStruct[0][0][iElm + 1, 0]
-                viscCgWing[iElm, 1] = self.viscStruct[0][0][iElm, 1] + self.viscStruct[0][0][iElm + 1, 1]
-                viscCgWing[iElm, 2] = self.viscStruct[0][0][iElm, 2] + self.viscStruct[0][0][iElm + 1, 2]
-            viscCgWing/=2
-
-            viscCgWake = np.zeros((len(self.viscStruct[0][1][:,0]) - 1, 3))
-            for iElm in range(len(self.viscStruct[0][1][:,0]) - 1):
-                viscCgWake[iElm, 0] = self.viscStruct[0][1][iElm, 0] + self.viscStruct[0][1][iElm + 1, 0]
-                viscCgWake[iElm, 1] = self.viscStruct[0][1][iElm, 1] + self.viscStruct[0][1][iElm + 1, 1]
-                viscCgWake[iElm, 2] = self.viscStruct[0][1][iElm, 2] + self.viscStruct[0][1][iElm + 1, 2]
-            viscCgWake/=2
-
-            invCgWing = np.zeros((len(self.bndWing.tag.elems), 3))
-            for iElm, e in enumerate(self.bndWing.tag.elems):
-                for n in e.nodes:
-                    invCgWing[iElm, 0] +=n.pos[0]
-                    invCgWing[iElm, 1] +=n.pos[1]
-                    invCgWing[iElm, 2] +=n.pos[2]
-                invCgWing[iElm, :] /= e.nodes.size()
-            invCgWake = np.zeros((len(self.bndWake.tag.elems), 3))
-            for iElm, e in enumerate(self.bndWake.tag.elems):
-                for n in e.nodes:
-                    invCgWake[iElm, 0] +=n.pos[0]
-                    invCgWake[iElm, 1] +=n.pos[1]
-                    invCgWake[iElm, 2] +=n.pos[2]
-                invCgWake[iElm, :] /= e.nodes.size()
-            
-            """plt.plot(self.invStruct[0][:,0], self.invStruct[0][:,1], 'x')
-            plt.plot(invCgWing[:, 0], invCgWing[:,1], 'o')
-            plt.show()"""
-
-        self.tms['Interpolation'].start()
-        mappedBlwWing = RBFInterpolator(viscCgWing, BlwWing, neighbors=5, kernel='linear', smoothing=1e-6)(invCgWing)
-        mappedBlwWake = RBFInterpolator(viscCgWake, BlwWake, neighbors=5, kernel='linear', smoothing=1e-6)(invCgWake)
-        self.tms['Interpolation'].stop()
-        for iElm in range(len(self.bndWing.tag.elems)):
-            self.bndWing.setU(iElm, mappedBlwWing[iElm])
-        for iElm in range(len(self.bndWake.tag.elems)):
-            self.bndWake.setU(iElm, mappedBlwWake[iElm])
-
-
-    def RunSolver(self):
-        self.iSolver.run()
-
-    def ResetSolver(self):
-        self.iSolver.reset()
-
-    def GetCp(self,bnd):
-        import dart.utils as floU
-
-        Cp = floU.extract(bnd, self.iSolver.Cp)
-        return Cp
-
-    def GetAlpha(self):
-        return self.iSolver.pbl.alpha
-
-    def GetMinf(self):
-        return self.iSolver.pbl.M_inf
-    
-    def GetCl(self):
-        return self.iSolver.Cl
-    
-    def GetCm(self):
-        return self.iSolver.Cm
-
-    def __Remesh(self, vec, idx_stag):
-        xUp = vec[:idx_stag+1]
-        xUp = xUp[::-1]
-
-        xLw = vec[idx_stag:]
-        return xUp, xLw
-
-    def __RbfToSections(self, var, varWk):
-        varWk = np.reshape(varWk, [len(varWk), 1])
-        completeWake = np.hstack((self.invStruct[1], varWk))
-        _, idx = np.unique(completeWake[:,:3], axis=0, return_index=True)
-        completeWake = completeWake[np.sort(idx)]
-
-        var = np.reshape(var, [len(var), 1])
-        completeWing = np.hstack((self.invStruct[0], var))
-
-        mappedVar = [[np.zeros(0) for _ in range(3)] for _ in range(len(self.viscStruct))]
-
-        for iSec in range(len(mappedVar)):
-            for iReg in range(len(mappedVar[iSec])):
-                if iReg==0: # Wing
-                    mappedVar[iSec][iReg] = RBFInterpolator(completeWing[:,:3], completeWing[:,3], neighbors=5, kernel='linear', smoothing=1e-2)(self.viscStruct[iSec][0])
-                if iReg==1: # Wake
-                    mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:3], completeWake[:,3], neighbors=5, kernel='linear', smoothing=1e-2)(self.viscStruct[iSec][1])
-        
-        """plt.plot(completeWing[:,0], completeWing[:,3], 'x')
-        plt.plot(completeWake[:,0], completeWake[:,3], 'x')
-        plt.plot(self.viscStruct[0][0][:,0], mappedVar[0][0])
-        plt.plot(self.viscStruct[0][1][:,0], mappedVar[0][1])
-        plt.show()"""
-        return mappedVar
-
-
-    def __GetWing(self, bnd_wing, bnd_wake, config):
-
-        if config['nDim'] == 2 and config['nSections'] != 1:
-            raise RuntimeError('Cannot create more than 1 section on 2D geometry. Specified:', config['nSections'])
-
-        viscStruct = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0], 3)) for i in range(2)] for _ in range(config['nSections'])]
-
-        invStruct = [np.zeros((len(bnd_wake.nodes) if iReg==1 else len(bnd_wing.nodes), 3)) for iReg in range(2)]
-
-        for iNode, n in enumerate(bnd_wing.nodes):
-            invStruct[0][iNode, 0] = n.pos[0]
-            invStruct[0][iNode, 1] = n.pos[1]
-            invStruct[0][iNode, 2] = n.pos[2]
-        for iNode, n in enumerate(bnd_wake.nodes):
-            invStruct[1][iNode, 0] = n.pos[0]
-            invStruct[1][iNode, 1] = n.pos[2]
-            invStruct[1][iNode, 2] = n.pos[1]
-
-        if config['nDim'] == 2:
-            spanDistri = np.zeros(1)
-            spanDistri[0] = invStruct[0][0,2] # Span direction is the z direction 
-
-        elif config['nDim'] == 3:
-            invStruct[0][:,[1,2]] = invStruct[0][:,[2,1]] # Inverse y and z so that spanwise direction is in z
-            invStruct[1][:,[1,2]] = invStruct[1][:,[2,1]] # Inverse y and z so that spanwise direction is in z
-            spanDistri = np.linspace(min(invStruct[0][:,2]), max(invStruct[0][:,2]), config['nSections'])
-            spanDistri[0] +=0.1
-            spanDistri[-1] -=0.1
-
-        for iSec, z in enumerate(spanDistri):
-
-            # Create x distribution at the corresponding span position
-            minX = min(invStruct[0][abs(invStruct[0][:,2]-z)<=config['eps'], 0])
-            chordLength = max(invStruct[0][abs(invStruct[0][:,2]-z)<=config['eps'], 0]) - minX
-            xDistri = minX + 1/2*(1 + np.cos(np.linspace(0, np.pi, config['nPoints'][0])))*chordLength
-            nElmSide = len(xDistri)
-            viscStruct[iSec][0][:,0] = np.concatenate((xDistri, xDistri[::-1]), axis=None)
-            viscStruct[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
-
-            wakeLength = np.max(invStruct[1][:,0]) - np.max(viscStruct[iSec][0][:,0])
-            viscStruct[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
-            viscStruct[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
-
-            if config['nDim'] == 2:
-                upper, lower = self.__SeparateArf(invStruct[0][:,0], invStruct[0][:,1])
-                viscStruct[iSec][0][:nElmSide,1] = interp1d(upper[:,0], upper[:,1])(viscStruct[iSec][0][:nElmSide,0])
-                viscStruct[iSec][0][nElmSide:,1] = interp1d(lower[:,0], lower[:,1])(viscStruct[iSec][0][nElmSide:,0])
-                viscStruct[iSec][0] = np.delete(viscStruct[iSec][0], nElmSide, axis=0)
-
-                viscStruct[iSec][1][:,1] = interp1d(invStruct[1][:,0], invStruct[1][:,1])(viscStruct[iSec][1][:,0])
-
-                """plt.plot(viscStruct[iSec][0][:,0], viscStruct[iSec][0][:,1], '-')
-                plt.plot(viscStruct[iSec][0][:nElmSide,0], viscStruct[iSec][0][:nElmSide,1], 'x')
-                plt.plot(viscStruct[iSec][0][nElmSide:,0], viscStruct[iSec][0][nElmSide:,1], 'o')
-                plt.plot(viscStruct[iSec][1][:,0], np.zeros(len(viscStruct[iSec][1][:,0])), 'x-')
-                plt.xlim([-0.5,2])
-                plt.ylim([-0.5, 0.5])
-                plt.show()"""
-
-            elif config['nDim'] == 3:
-                nPointsMin = (config['k'][0]+1) * (config['k'][1]+1) # See Scipy bisplrep doc
-
-                wingSpan = max(invStruct[0][:,2]) - min(invStruct[0][:,2])
-                portion = invStruct[0][abs(invStruct[0][:,2] - z) < 0.01 * wingSpan,:]
-
-                portionUpper, portionLower = self.__SeparateArf(portion[:,0], portion[:,1], portion[:,2])
-                if len(portionUpper)<=nPointsMin or len(portionLower)<=nPointsMin:
-                    print(ccolors.ANSI_YELLOW,'dart::Interpolator Warning: Found {:.0f} points within 1% of spanwise position z={:.2f}. Trying 10%. '.format(
-                        len(portion), z),ccolors.ANSI_RESET)
-                    portion = invStruct[0][abs(invStruct[0][:,2] - z) < 0.1 * wingSpan,:]
-                    if len(portion)<=nPointsMin:
-                        print(ccolors.ANSI_RED, 'Failed. Continuing\n', ccolors.ANSI_RESET)
-                    else:
-                        print(ccolors.ANSI_GREEN, 'Sucess.\n', ccolors.ANSI_RESET)
-
-                # Upper
-                sWing = (len(portionUpper[:,0])-np.sqrt(2*len(portionUpper[:,0])) + len(portionUpper[:,0])+np.sqrt(2*len(portionUpper[:,0])))/2
-                tckWingUp = bisplrep(portionUpper[:,0], portionUpper[:,2], portionUpper[:,1], kx=config['k'][0], ky=config['k'][1], s=sWing)
-                for iPoint in range(nElmSide):
-                    viscStruct[iSec][0][iPoint,1] = bisplev(viscStruct[iSec][0][iPoint,0], viscStruct[iSec][0][iPoint,2], tckWingUp)
-                
-                # Lower
-                sWing = (len(portionLower[:,0])-np.sqrt(2*len(portionLower[:,0])) + len(portionLower[:,0])+np.sqrt(2*len(portionLower[:,0])))/2
-                tckWingUp = bisplrep(portionLower[:,0], portionLower[:,2], portionLower[:,1], kx=config['k'][0], ky=config['k'][1], s=sWing)
-                for iPoint in range(nElmSide, len(viscStruct[iSec][0][:,0])):
-                    viscStruct[iSec][0][iPoint,1] = bisplev(viscStruct[iSec][0][iPoint,0], viscStruct[iSec][0][iPoint,2], tckWingUp)
-
-                # Wake
-                sWake = (len(invStruct[1][:,0])-np.sqrt(2*len(invStruct[1][:,0])) + len(invStruct[1][:,0])+np.sqrt(2*len(invStruct[1][:,0])))/2
-                tckWingUp = bisplrep(invStruct[1][:,0], invStruct[1][:,2], invStruct[1][:,1], kx=config['k'][0], ky=config['k'][1], s=sWake)
-                for iPoint in range(len(viscStruct[iSec][1][:,0])):
-                    viscStruct[iSec][1][iPoint,1] = bisplev(viscStruct[iSec][1][iPoint,0], viscStruct[iSec][1][iPoint,2], tckWingUp)
-
-
-        if config['nDim'] ==3:
-            fig = plt.figure()
-            ax = fig.add_subplot(projection='3d')
-            ax.scatter(invStruct[0][:,0], invStruct[0][:,1], invStruct[0][:,2], s=0.2)
-            for iSec in range(len(viscStruct)):
-                ax.scatter(viscStruct[iSec][0][:,0], viscStruct[iSec][0][:,1], viscStruct[iSec][0][:,2], s=0.9, color='red')
-            ax.set_ylim([-0.5, 0.5])
-            ax.set_xlabel('x')
-            ax.set_ylabel('y')
-            ax.set_zlabel('z')
-            plt.show()
-        return invStruct, viscStruct
-
-    def __SeparateArf(self, xInp, yInp, zInp=None):
-        if zInp is None:
-            zInp = np.zeros(len(xInp))
-        if len(xInp) != len(yInp) or len(xInp) != len(zInp):
-            raise RuntimeError ('dart::Interpolator::__SeparateArf Missmatch between intput vector sizes.')
-
-        Q=np.zeros((len(xInp), 3))
-        Q[:,0] = xInp.copy()
-        Q[:,1] = yInp.copy()
-        Q[:,2] = zInp.copy()
-
-        Q = Q[Q[:, 0].argsort()]
-
-        le = Q[0,:]
-        te = Q[-1,:]
-        ymean = (le[1] + te[1]) / 2
-
-        upper, lower, save = le, le, le
-
-        flag1 = 0
-        iPoint = 0
-        for i in range(len(Q[:,0])):
-            if Q[iPoint,1] == 0.0:
-                if flag1 == 1:
-                    Q=np.delete(Q, iPoint, axis = 0)
-                    iPoint-=1
-                else:
-                    flag1 = 1
-            iPoint +=1
-
-        for iPoint in range(1, len(Q[:,0])):
-            if Q[iPoint,1] > ymean:
-                upper = np.vstack((upper, Q[iPoint,:]))
-            elif Q[iPoint, 1] < ymean:
-                lower = np.vstack((lower, Q[iPoint,:]))
-            else:
-                save = np.vstack((save, Q[iPoint, :]))
-        
-        _, idx = np.unique(save, return_index=True, axis=0)
-        save=save[np.sort(idx)]
-
-        #save = np.delete(save, 0, axis=0)
-
-        
-        try:
-            if len(save[:,0])>0:
-                upper = np.vstack((upper, save))
-                lower = np.vstack((lower, save))
-        except:
-            pass
-
-        upper = np.vstack((upper, te))
-        lower = np.vstack((lower, te))
-
-        return upper, lower
\ No newline at end of file
diff --git a/dart/pyVII/vRbf.py b/dart/pyVII/vRbf.py
new file mode 100644
index 0000000000000000000000000000000000000000..d227ef759379b44f90d776e5cc40f662d586dc0c
--- /dev/null
+++ b/dart/pyVII/vRbf.py
@@ -0,0 +1,134 @@
+from audioop import lin2adpcm
+import numpy as np
+
+class RbfInterpolator:
+
+    def __init__(self) -> None:
+        pass
+
+    def RBFInterp(self, xs, ys, x, type, R):
+        """"
+        xs : numpy.ndarray Nx3, base mesh
+        ys : numpy.ndarray Nx1, solution (on xs)
+        x : numpy.ndarray Mx3, target mesh
+        type : string
+        R : double Radius"""
+        if len(ys) != len(xs[:,0]):
+            raise RuntimeError ("vector size does not match mesh", len(ys), "expected:",len(xs[:,0]))
+
+        fPar = self.RBFparam(xs, ys, type, R)
+
+        #if ~isempty(fPar)
+        # Evaluate interpolation function
+        y = self.RBFeval(xs, x, fPar, type, R);
+        #else
+        #    y = [];
+        return y
+
+    def RBFParam(self, xs, ys, RBFtype, R):
+
+        if len(xs[:, 0]) == len(ys[:,0]):
+            
+            Ns = len(xs[:,0])
+            Ms = len(ys)
+            dim = len(xs[0,:])
+
+            ncol = 1 + dim
+            r = np.zeros(Ns)
+
+            for i in range(Ns):
+                for j in range(i,Ns):
+                    r[i, j] = np.norm( xs[i, :] - xs[j, :])
+                    r[j, i] = r[i, j]
+
+            P = [np.ones((Ns, 1)), xs]
+            M = self.radialFunction(r, RBFtype, R)
+            
+            A = [[M, P], [np.transpose(P), np.zeros(ncol)]]
+            b = [[ys], [np.zeros(ncol, Ms)]]
+
+            fPar = np.linalg.solve(A,b)
+            
+        else:
+            
+            fPar = []
+            M = []
+        return fPar
+
+    def radialFunction(self, r, RBFtype, R):
+
+        r = r/R;
+
+        phi = np.zeros(len(r))
+
+        if RBFtype == 'R1':
+            phi = r
+        elif RBFtype == 'R3':
+            phi = r**3
+        elif RBFtype == 'TPS2':
+            phi[r>0] = r[r>0]**2*np.log(r[r>0])
+        elif RBFtype == 'Q':
+            phi = 1 + r**2
+        elif RBFtype == 'MQ':
+            phi = np.sqrt(1 + r**2)
+        elif RBFtype == 'IMQ':
+            phi = 1/np.sqrt(1 + r**2)
+        elif RBFtype == 'IQ':
+            phi = 1/(1 + r**2)
+        elif RBFtype == 'GS':
+            phi = np.exp(-r**2)
+        elif RBFtype == 'CP_C0':
+            I = (r < 1);
+            phi[r<1] = (1 - r[r<1])**2
+        elif RBFtype == 'CP_C2':
+            I = (r < 1)
+            phi[r<1] = (1 - r[r<1])**4*(4*r[r<1]+ 1)
+        elif RBFtype == 'CP_C4':
+            I = (r < 1)
+            phi[r<1] = (1 - r(I))**6*(35/3*r[r<1]**2 + 6*r[r<1] + 1)
+        """if RBFtype == 'CP_C6':
+            I = (r < 1)
+            phi(I) = (1 - r(I)).^8.*(32*r(I).^3 + 25*r(I).^2 + 8*r(I) + 1)
+        if RBFtype == 'CTPS_C0':
+            I = (r < 1)
+            phi(I) = (1 - r(I)).^5
+        if RBFtype == 'CTPS_C1':
+            I = (r < 1 & r > 0)
+            phi(I) = 1 + 80/3*r(I).^2 - 40*r(I).^3 + 15*r(I).^4 - 8/3*r(I).^5 + 20*r(I).^2.*log(r(I))
+            phi(r == 0) = 1;
+        if RBFtype == 'CTPS_C2a':
+            I = (r < 1 & r > 0)
+            phi(I) = 1 - 30*r(I).^2 - 10*r(I).^3 + 45*r(I).^4 - 6*r(I).^5 - 60*r(I).^3.*log(r(I))
+            phi(r == 0) = 1
+        if RBFtype == 'CTPS_C2b':
+            I = (r < 1 & r > 0)
+            phi(I) = 1 - 20*r(I).^2 + 80*r(I).^3 - 45*r(I).^4 -16*r(I).^5 + 60*r(I).^4.*log(r(I))
+            phi(r == 0) = 1
+        else:
+            phi = r"""
+
+        return phi
+
+    def RBFeval(self, xs, x, fPar, RBFtype, R):
+
+        Ns = len(xs[:,0]);
+        dim = len(xs[0,:]);
+
+        if len(xs[0,:]) == dim and len(fPar[0,:]) == Ns + dim + 1:
+            
+            N = len(x[:,0])
+            
+            r = np.zeros((N, Ns))
+            for i in range(N):
+                for j in range(Ns):
+                    r[i, j] = np.norm( x[i, :] - xs[j, :] )
+            
+            P = [np.ones((N, 1)), x]
+            M = self.radialFunction(r, RBFtype, R)
+            
+            y = [M, P]*fPar;
+            
+        else:
+            
+            y = []
+        return y
\ No newline at end of file
diff --git a/dart/src/wAdjoint.cpp b/dart/src/wAdjoint.cpp
index 00e62aa3463dcd1ffe109e2bd7c81b572e91399d..c3a534016cec1e2b928bed802151dd13a617f5a5 100644
--- a/dart/src/wAdjoint.cpp
+++ b/dart/src/wAdjoint.cpp
@@ -662,4 +662,4 @@ void Adjoint::write(std::ostream &out) const
 {
     out << "dart::Adjoint"
         << "\n";
-}
+}
\ No newline at end of file
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index 5cf782b8d8d2bae78bcc6409e7c5844903cbb631..c1ed9ac5717d3cb311712ce3fb061a8eb8917baf 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -36,7 +36,6 @@ void BLRegion::SetMesh(std::vector<double> x,
 
     if (mesh->GetDiffnElm()!=0)
     {
-        std::cout << "Resizing vectors" << std::endl;
         U.resize(nVar * nMarkers, 0);
         Ret.resize(nMarkers, 0);
         Cd.resize(nMarkers);
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 84f5388d6e5d7bc8a42e77456a3bb991317f9eb0..6762667511945d2018035a55e4c49afc18d51721 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -190,7 +190,6 @@ int ViscSolver::Run(unsigned int couplIter)
                 std::vector<double> LowerCond(Sections[iSec][iRegion]->GetnVar(), 0);
                 UpperCond = std::vector<double>(Sections[iSec][0]->U.end() - Sections[iSec][0]->GetnVar(), Sections[iSec][0]->U.end()); /* Base std::vector constructor to slice */
                 LowerCond = std::vector<double>(Sections[iSec][1]->U.end() - Sections[iSec][1]->GetnVar(), Sections[iSec][1]->U.end()); /* Base std::vector constructor to slice */
-                std::cout << " Size uppercond " << UpperCond.size()  << " Size lowercond " << LowerCond.size() << std::endl;
                 Sections[iSec][iRegion]->SetWakeBC(Re, UpperCond, LowerCond);
                 lockTrans = true;
             }
@@ -204,7 +203,7 @@ int ViscSolver::Run(unsigned int couplIter)
                 tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
 
                 /* Time integration */
-                if (iRegion == 1 && iPoint == 1)
+                /*if (iRegion == 1 && iPoint == 1)
                 {
                     std::scientific;
                 std::cout << " iSec " << iSec
@@ -216,15 +215,10 @@ int ViscSolver::Run(unsigned int couplIter)
                           << " DeltaStarExt " << Sections[iSec][iRegion]->invBnd->GetDeltaStar(iPoint)
                           << std::endl;
                 Sections[iSec][iRegion]->PrintSolution(iPoint);
-                }
+                }*/
 
                 pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
 
-                if (Sections[iSec][iRegion]->Cf[iPoint] < 0)
-                {
-                    std::cout << "Negative Cf at point " << iPoint << " : " << Sections[iSec][iRegion]->Cf[iPoint] << std::endl;
-                }
-
                 /* if (iRegion == 1 && iPoint == 265)
                 {
                 std::cout << " iSec " << iSec
diff --git a/dart/tests/bli3.py b/dart/tests/bli3.py
index 8083eee7f8e359e237658011d7a5e08e70e7a417..2f0a77e38711fc95260b22cb377da9a543a6ea51 100644
--- a/dart/tests/bli3.py
+++ b/dart/tests/bli3.py
@@ -31,7 +31,7 @@ import dart.default as floD
 
 import dart.pyVII.vUtils as viscU
 import dart.pyVII.vCoupler as viscC
-import dart.pyVII.vInterpolator2 as vInterpol
+import dart.pyVII.vInterpolator as vInterpol
 
 import fwk
 from fwk.testing import *
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index eddabfa0ad56bdb97d12f89d00dda8498c442ff6..ce4f843ff16942ac4a8d8a0d6c0db878b580dc15 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python3 
 # -*- coding: utf-8 -*-
 
 # Copyright 2020 University of Liège
@@ -39,6 +39,8 @@ import dart.default as floD
 #import dart.viscUtils as viscU
 
 import dart.pyVII.vCoupler as floVC
+import dart.pyVII.vCoupler2 as viscC # Rbf coupler
+import dart.pyVII.vInterpolator as Interpol # Rbf interpolator
 import numpy as np
 
 import tbox
@@ -57,9 +59,9 @@ def main():
 
     # define flow variables
     Re = 1e7
-    alpha = 2*math.pi/180
+    alpha = 5*math.pi/180
     #alphaSeq = np.arange(-5, 5.5, 0.5)
-    M_inf = 0.715
+    M_inf = 0.
 
     meshFactor = 1
     CFL0 = 1
@@ -67,6 +69,7 @@ def main():
     plotVar = 'H'
 
     nSections = 1
+    rbf = 0
 
     # ========================================================================================
 
@@ -76,8 +79,8 @@ def main():
     # mesh the geometry
     print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
     tms['msh'].start()
-    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.005}
-    msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
+    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
     # set the problem and add medium, initial/boundary conditions
@@ -93,9 +96,14 @@ def main():
     isolver = floD.newton(pbl)
     vsolver = floD.initBL(Re, M_inf, CFL0, nSections)
 
-    coupler = floVC.Coupler(isolver, vsolver)
 
-    coupler.CreateProblem(blw[0], blw[1])
+    if rbf:
+        config = {'nDim': dim, 'nSections':nSections, 'nPoints':[100, 50], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 10}
+        iSolverAPI = Interpol.Interpolator(isolver, blw[0], blw[1], config)
+        coupler = viscC.Coupler(iSolverAPI, vsolver)
+    else:
+        coupler = floVC.Coupler(isolver, vsolver)
+        coupler.CreateProblem(blw[0], blw[1])
 
     #coupler.RunPolar(alphaSeq)
     coupler.Run()