diff --git a/CMakeCache.txt b/CMakeCache.txt
new file mode 100644
index 0000000000000000000000000000000000000000..dc789f0516728e350d38ef799a00b75f41ca59e9
--- /dev/null
+++ b/CMakeCache.txt
@@ -0,0 +1,38 @@
+# This is the CMakeCache file.
+# For build in directory: /Users/paulzer/lab/dartflo
+# It was generated by CMake: /usr/local/Cellar/cmake/3.22.3/bin/cmake
+# You can edit this file to change values found and used by cmake.
+# If you do not want to change any of the values, simply exit the editor.
+# If you do want to change a value, simply edit, save, and exit the editor.
+# The syntax for the file is as follows:
+# KEY:TYPE=VALUE
+# KEY is the name of a variable in the cache.
+# TYPE is a hint to GUIs for the type of VALUE, DO NOT EDIT TYPE!.
+# VALUE is the current value for the KEY.
+
+########################
+# EXTERNAL cache entries
+########################
+
+
+########################
+# INTERNAL cache entries
+########################
+
+//This is the directory where this CMakeCache.txt was created
+CMAKE_CACHEFILE_DIR:INTERNAL=/Users/paulzer/lab/dartflo
+//Major version of cmake used to create the current loaded cache
+CMAKE_CACHE_MAJOR_VERSION:INTERNAL=3
+//Minor version of cmake used to create the current loaded cache
+CMAKE_CACHE_MINOR_VERSION:INTERNAL=22
+//Patch version of cmake used to create the current loaded cache
+CMAKE_CACHE_PATCH_VERSION:INTERNAL=3
+//Path to CMake executable.
+CMAKE_COMMAND:INTERNAL=/usr/local/Cellar/cmake/3.22.3/bin/cmake
+//Path to cpack program executable.
+CMAKE_CPACK_COMMAND:INTERNAL=/usr/local/Cellar/cmake/3.22.3/bin/cpack
+//Path to ctest program executable.
+CMAKE_CTEST_COMMAND:INTERNAL=/usr/local/Cellar/cmake/3.22.3/bin/ctest
+//Path to CMake installation.
+CMAKE_ROOT:INTERNAL=/usr/local/Cellar/cmake/3.22.3/share/cmake
+
diff --git a/CMakeFiles/cmake.check_cache b/CMakeFiles/cmake.check_cache
new file mode 100644
index 0000000000000000000000000000000000000000..3dccd731726d7faa8b29d8d7dba3b981a53ca497
--- /dev/null
+++ b/CMakeFiles/cmake.check_cache
@@ -0,0 +1 @@
+# This file is generated by cmake for dependency checking of the CMakeCache.txt file
diff --git a/dart/default.py b/dart/default.py
index 3a03df5626170c58ef3be646579d518954072a69..b697a3f92f7ec2386c0d87b922212db0305df4f8 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -23,6 +23,11 @@ from fwk.wutils import parseargs
 import dart
 import tbox
 
+def initBL(Re, Minf, CFL0, meshFactor):
+    solver = dart.ViscSolver(Re, Minf, CFL0, meshFactor)
+    return solver
+
+
 def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
     """Initialize mesh and mesh writer
 
diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index 79f4459069fb17b129aba63c8c5a5af758a87e3e..449f2a2ce2c669f6c4e5744220c3d2a5eae57e03 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -19,54 +19,61 @@
 ## VII Coupler
 # Paul Dechamps
 
+from dart.viscous.Master.DAirfoil import Airfoil
+from dart.viscous.Master.DWake import Wake
 from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
 
 from fwk.coloring import ccolors
 import math
 
 import fwk
-import dart.pyVII.vUtils as viscU
+#import dart.viscUtils 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
+  def __init__(self, iSolver, vSolver) -> None:
+    self.iSolver=iSolver
+    self.vSolver=vSolver
 
     self.maxCouplIter = 150
     self.couplTol = 1e-4
 
     self.tms=fwk.Timers()
-    self.resetInv = True
+    self.resetInv = False
+    pass
+
+  def CreateProblem(self, _boundaryAirfoil, _boundaryWake):
+
+    self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects
+    self.blwVel = [None, None, None]
+    self.deltaStar = [None, None, None]
+    self.xx = [None, None, None]
     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('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'
+    .format(self.iSolver.pbl.alpha*180/math.pi, self.iSolver.pbl.M_inf, self.vSolver.Re/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.iSolver.reset()
+      self.iSolver.run()
       self.tms['inviscid'].stop()
 
       # Extract inviscid boundary
       self.tms['processing'].start()
-      self.iSolverAPI.GetInvBC(self.vSolver, couplIter)
+      self.__ImposeInvBC(couplIter)
       self.tms['processing'].stop()
 
       # Run viscous solver
@@ -76,35 +83,32 @@ class Coupler:
       # print('Viscous ops terminated with exit code ', exitCode)
 
       self.tms['processing'].start()
-      self.iSolverAPI.SetBlowingVel(self.vSolver)
+      self.__ImposeBlwVel()
       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:
+      if error  <= self.couplTol and exitCode==0:
 
-        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
+        print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>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
 
 
-      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)))
+      if couplIter%5==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,
+        self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[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
+      return couplIter
 
 
   def RunPolar(self, alphaSeq):
@@ -127,4 +131,102 @@ class Coupler:
       coeffTable[1,ind] = self.iSolver.Cl
       coeffTable[2,ind] = self.vSolver.Cdt
 
-    return coeffTable
\ No newline at end of file
+    return coeffTable
+
+  def __ImposeBlwVel(self):
+
+    # Retreive and set blowing velocities.
+
+    for iRegion in range(3):
+      self.blwVel[iRegion] = np.asarray(self.vSolver.ExtractBlowingVel(0, iRegion))
+      self.deltaStar[iRegion] = np.asarray(self.vSolver.ExtractDeltaStar(0, iRegion))
+      self.xx[iRegion] = np.asarray(self.vSolver.Extractxx(0, iRegion))
+
+    blwVelAF = np.concatenate((self.blwVel[0], self.blwVel[1]))
+    deltaStarAF = np.concatenate((self.deltaStar[0], self.deltaStar[1][1:])) # Remove stagnation point doublon.
+    xxAF = np.concatenate((self.xx[0], self.xx[1][1:])) # Remove stagnation point doublon.
+    blwVelWK = self.blwVel[2]
+    deltaStarWK = self.deltaStar[2]
+    xxWK = self.xx[2]
+
+    self.group[0].u = blwVelAF[self.group[0].connectListElems.argsort()]
+    self.group[1].u = blwVelWK[self.group[1].connectListElems.argsort()]
+
+    self.group[0].deltaStar = deltaStarAF[self.group[0].connectListNodes.argsort()]
+    self.group[1].deltaStar = deltaStarWK[self.group[1].connectListNodes.argsort()]
+
+    self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()]
+    self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()]
+
+    for n in range(0, len(self.group)):
+      for i in range(0, self.group[n].nE):
+        self.group[n].boundary.setU(i, self.group[n].u[i])
+
+  def __ImposeInvBC(self, couplIter):
+
+    for n in range(0, len(self.group)):
+
+      for i in range(0, len(self.group[n].boundary.nodes)):
+
+        self.group[n].v[i,0] = self.iSolver.U[self.group[n].boundary.nodes[i].row][0]
+        self.group[n].v[i,1] = self.iSolver.U[self.group[n].boundary.nodes[i].row][1]
+        self.group[n].v[i,2] = self.iSolver.U[self.group[n].boundary.nodes[i].row][2]
+        self.group[n].Me[i] = self.iSolver.M[self.group[n].boundary.nodes[i].row]
+        self.group[n].rhoe[i] = self.iSolver.rho[self.group[n].boundary.nodes[i].row]
+    
+    dataIn=[None,None]
+    for n in range(0,len(self.group)):
+      self.group[n].connectListNodes, self.group[n].connectListElems, dataIn[n]  = self.group[n].connectList()
+    fMeshDict, _, _ = GroupConstructor().mergeVectors(dataIn, 0, 1)
+
+    if ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) != self.nElmUpperPrev) or couplIter==0):
+
+      # print('STAG POINT MOVED')
+      
+      # Initialize mesh
+      self.vSolver.SetGlobMesh(0, 0, fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['yGlobCoord'][:fMeshDict['bNodes'][1]], fMeshDict['zGlobCoord'][:fMeshDict['bNodes'][1]])
+      self.vSolver.SetGlobMesh(0, 1, fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['yGlobCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['zGlobCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+      self.vSolver.SetGlobMesh(0, 2, fMeshDict['globCoord'][fMeshDict['bNodes'][2]:], fMeshDict['yGlobCoord'][fMeshDict['bNodes'][2]:], fMeshDict['zGlobCoord'][fMeshDict['bNodes'][2]:])
+      #self.vSolver.InitMeshes(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['globCoord'][:fMeshDict['bNodes'][1]])
+      #self.vSolver.InitMeshes(1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+      #self.vSolver.InitMeshes(2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:], fMeshDict['globCoord'][fMeshDict['bNodes'][2]:])
+
+      self.vSolver.SetxxExt(0, 0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]])
+      self.vSolver.SetxxExt(0, 1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])    
+      self.vSolver.SetxxExt(0, 2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:])
+      
+      self.vSolver.SetDeltaStarExt(0, 0, fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
+      self.vSolver.SetDeltaStarExt(0, 1, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])    
+      self.vSolver.SetDeltaStarExt(0, 2, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
+
+    else:
+
+      self.vSolver.SetxxExt(0, 0, fMeshDict['xxExt'][:fMeshDict['bNodes'][1]])
+      self.vSolver.SetxxExt(0, 1, fMeshDict['xxExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])    
+      self.vSolver.SetxxExt(0, 2, fMeshDict['xxExt'][fMeshDict['bNodes'][2]:])
+      
+      self.vSolver.SetDeltaStarExt(0, 0, fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
+      self.vSolver.SetDeltaStarExt(0, 1, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])    
+      self.vSolver.SetDeltaStarExt(0, 2, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
+
+    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]])
+    self.vSolver.SetVelocities(0, 2, fMeshDict['vx'][fMeshDict['bNodes'][2]:], fMeshDict['vy'][fMeshDict['bNodes'][2]:], fMeshDict['vz'][fMeshDict['bNodes'][2]:])
+    #self.vSolver.SendInvBC(0, fMeshDict['vtExt'][:fMeshDict['bNodes'][1]], fMeshDict['Me'][:fMeshDict['bNodes'][1]], fMeshDict['rho'][:fMeshDict['bNodes'][1]])
+    #self.vSolver.SendInvBC(1, fMeshDict['vtExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+    #self.vSolver.SendInvBC(2, fMeshDict['vtExt'][fMeshDict['bNodes'][2]:], fMeshDict['Me'][fMeshDict['bNodes'][2]:], fMeshDict['rho'][fMeshDict['bNodes'][2]:])
+    self.vSolver.SetMe(0, 0, fMeshDict['Me'][:fMeshDict['bNodes'][1]])
+    self.vSolver.SetMe(0, 1, fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+    self.vSolver.SetMe(0, 2, fMeshDict['Me'][fMeshDict['bNodes'][2]:])
+    self.vSolver.SetRhoe(0, 0, fMeshDict['rho'][:fMeshDict['bNodes'][1]])
+    self.vSolver.SetRhoe(0, 1, fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+    self.vSolver.SetRhoe(0, 2, fMeshDict['rho'][fMeshDict['bNodes'][2]:])
+    # print('Inviscid boundary imposed')
+  
+
+
+  
+
+
+
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 6bc2ac7dfb3f26daac90a782af0b12784c5c4532..586cdacfd4e989d1ff2f9c82a660464871ca8cf0 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -433,7 +433,7 @@ void ViscSolver::SetGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x
 
 void ViscSolver::SetVelocities(size_t iSec, size_t iRegion, std::vector<double> vx, std::vector<double> vy, std::vector<double> vz)
 {
-    std::vector<double> Vt(vx.size());
+    std::vector<double> Vt(vx.size(), 0.0);
     double alpha=0;
 
     for (size_t iPoint=0; iPoint<vx.size()-1; ++iPoint)
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index d827250aafd2650a4a28dd2cce477b95c2811ba5..78b802f30cc977d041c3eeb115481fb5b8ead3bf 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -34,14 +34,11 @@
 # This test is provided to ensure that the solver works properly.
 # Mesh refinement may have to be performed to obtain physical results.
 
-from tkinter import NS
+import dart.utils as floU
 import dart.default as floD
+#import dart.viscUtils as viscU
 
-import dart.pyVII.vUtils as viscU
-import dart.pyVII.vCoupler as viscC
-import dart.pyVII.vInterpolator2 as vInterpol
-
-from matplotlib import pyplot as plt
+import dart.pyVII.vCoupler as floVC
 import numpy as np
 
 import tbox
@@ -60,19 +57,22 @@ def main():
 
     # define flow variables
     Re = 1e7
-    alpha = 12.*math.pi/180
+    alpha = 12*math.pi/180
     #alphaSeq = np.arange(-5, 5.5, 0.5)
     M_inf = 0.
 
+    meshFactor = 1
     CFL0 = 1
 
+    plotVar = 'H'
+
+    nSections = 1
+
     # ========================================================================================
 
     c_ref = 1
     dim = 2
-    
-    nSections = 1
-        
+
     # mesh the geometry
     print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
     tms['msh'].start()
@@ -89,31 +89,25 @@ def main():
 
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
     tms['solver'].start()
+    
+    isolver = floD.newton(pbl)
+    vsolver = floD.initBL(Re, M_inf, CFL0, nSections)
 
-    config = {'nDim': dim, 'nSections':nSections, 'nPoints':[500, 50], 'eps':0.02}
+    coupler = floVC.Coupler(isolver, vsolver)
 
-    #iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw[0], blw[1], nSections, dim)
-    dragIter = [np.zeros(0) for i in range(1)]
-    for i in range(len(dragIter)):
-        iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw[0], blw[1], config)
-        vSolver = viscU.initBL(Re, M_inf, CFL0, nSections)
-        coupler = viscC.Coupler(iSolverAPI, vSolver)
+    coupler.CreateProblem(blw[0], blw[1])
 
     #coupler.RunPolar(alphaSeq)
-        dragIter[i] = coupler.Run()
-        plt.plot(dragIter[i])
-    plt.title('{:.0f} real. {:.2f} deg AoA'.format(len(dragIter), alpha*180/np.pi))
-    plt.xlabel('Iters')
-    plt.ylabel('$c_d$')
+    coupler.Run()
     tms['solver'].stop()
 
     # extract Cp
-    Cp = iSolverAPI.GetCp(bnd.groups[0].tag.elems)
+    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
     tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
     # 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(Re/1e6, M_inf, alpha*180/math.pi, iSolverAPI.GetCl(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.GetCm()))
+    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(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
 
     # display timers
     tms['total'].stop()
@@ -122,15 +116,14 @@ def main():
     print(tms)
     print('SOLVERS statistics')
     print(coupler.tms)
-    print(iSolverAPI.tms)
-    plt.show()
-
 
+    xtr=vsolver.Getxtr()
+    
     """# visualize solution and plot results
     #floD.initViewer(pbl)
-    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(iSolver.Cl, vSolver.Cdt, iSolver.Cm, 4), True)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True)
 
-    blScalars, blVectors = viscU.ExtractVars(vSolver)
+    blScalars, blVectors = viscU.ExtractVars(vsolver)
     wData = viscU.Dictionary(blScalars, blVectors, xtr)
     viscU.PlotVars(wData, plotVar)"""
 
@@ -139,19 +132,19 @@ def main():
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
     if Re == 1e7 and M_inf == 0. and alpha == 2.*math.pi/180:
-        tests.add(CTest('Cl', iSolverAPI.GetCl(), 0.2271, 5e-3)) # Xfoil 0.2271
-        tests.add(CTest('Cd', vSolver.Cdt, 0.00531, 1e-3, forceabs=True))
-        tests.add(CTest('Cdp', vSolver.Cdp, 0.0009, 1e-3, forceabs=True))
-        tests.add(CTest('xtr_top', vSolver.Getxtr(0,0), 0.20, 0.03, forceabs=True))
-        tests.add(CTest('xtr_bot', vSolver.Getxtr(0,1), 0.50, 0.03, forceabs=True))
+        tests.add(CTest('Cl', isolver.Cl, 0.2208, 5e-3))
+        tests.add(CTest('Cd', vsolver.Cdt, 0.00531, 1e-3, forceabs=True))
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0009, 1e-3, forceabs=True))
+        tests.add(CTest('xtr_top', xtr[0], 0.20, 0.03, forceabs=True))
+        tests.add(CTest('xtr_bot', xtr[1], 0.50, 0.03, forceabs=True))
     if Re == 1e7 and M_inf == 0. and alpha == 5.*math.pi/180:
-        tests.add(CTest('Cl', iSolverAPI.GetCl(), 0.5488, 5e-3))
-        tests.add(CTest('Cd', vSolver.Cdt, 0.0062, 1e-3, forceabs=True))
-        tests.add(CTest('Cdp', vSolver.Cdp,0.0018, 1e-3, forceabs=True))
-        tests.add(CTest('xtr_top', vSolver.Getxtr(0,0), 0.06, 0.03, forceabs=True))
-        tests.add(CTest('xtr_bot', vSolver.Getxtr(0,1), 0.74, 0.03, forceabs=True))
-    """else:
-        raise Exception('Test not defined for this flow')"""
+        tests.add(CTest('Cl', isolver.Cl, 0.5488, 5e-3))
+        tests.add(CTest('Cd', vsolver.Cdt, 0.0062, 1e-3, forceabs=True))
+        tests.add(CTest('Cdp', vsolver.Cdp,0.0018, 1e-3, forceabs=True))
+        tests.add(CTest('xtr_top', xtr[0], 0.06, 0.03, forceabs=True))
+        tests.add(CTest('xtr_bot', xtr[1], 0.74, 0.03, forceabs=True))
+    else:
+        raise Exception('Test not defined for this flow')
 
     tests.run()
 
diff --git a/dart/viscous/Drivers/DGroupConstructor.py b/dart/viscous/Drivers/DGroupConstructor.py
index fab1f757d6968f25339b20eba95e6f76dba6dcda..2d91256ab083469eb10fa3cff1bf7fbe6b2a09a6 100755
--- a/dart/viscous/Drivers/DGroupConstructor.py
+++ b/dart/viscous/Drivers/DGroupConstructor.py
@@ -97,11 +97,17 @@ class GroupConstructor():
         else:
 
             Mesh         = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][:,0], dataIn[1][:,0]))
+            Meshy         = np.concatenate((dataIn[0][0][:,1], dataIn[0][1][:,1], dataIn[1][:,1]))
+            Meshz         = np.concatenate((dataIn[0][0][:,2], dataIn[0][1][:,2], dataIn[1][:,2]))
             MeshbNodes   = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][:,0]) + len(dataIn[1][:,0])]
 
             xx           = np.concatenate((xx_up, xx_lw, xx_wk))
             vtExt        = np.concatenate((vtExt_up, vtExt_lw, vtExt_wk))
 
+            vx = np.concatenate((dataIn[0][0][:,3], dataIn[0][1][:,3], dataIn[1][:,3]))
+            vy = np.concatenate((dataIn[0][0][:,4], dataIn[0][1][:,4], dataIn[1][:,4]))
+            vz = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][:,5], dataIn[1][:,5]))
+
             alpha        = np.concatenate((alpha_up, alpha_lw, alpha_wk))
             dv           = np.concatenate((dv_up, dv_lw, dv_wk))
 
@@ -110,8 +116,8 @@ class GroupConstructor():
             dStarext     = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][:,7], dataIn[1][:,7]))
             xxExt        = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][:,8], dataIn[1][:,8]))
 
-            fMeshDict    = {'globCoord': Mesh, 'bNodes': MeshbNodes, 'locCoord': xx,
-                            'vtExt': vtExt, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'xxExt': xxExt, 'dv': dv}
+            fMeshDict    = {'globCoord': Mesh, 'yGlobCoord': Meshy, 'zGlobCoord': Meshz, 'bNodes': MeshbNodes, 'locCoord': xx,
+                            'vtExt': vtExt, 'vx': vx, 'vy': vy, 'vz': vz, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'xxExt': xxExt, 'dv': dv}
 
             cMeshDict    = fMeshDict.copy()