From 85dfd5e2bef4d86ded74f3abd49f6666c69bde21 Mon Sep 17 00:00:00 2001
From: Paul <paul.dechamps@student.uliege.be>
Date: Sat, 2 Jul 2022 23:17:06 +0200
Subject: [PATCH] The great departure of vii :(

---
 dart/_src/dartw.h                             |   8 -
 dart/_src/dartw.i                             |   9 -
 dart/src/dart.h                               |  10 -
 dart/src/wInterpolator.cpp                    |  18 --
 dart/src/wInterpolator.h                      |  22 --
 dart/tests/bli2.py                            | 158 --------------
 dart/tests/bli_Polar.py                       | 124 -----------
 dart/tests/bli_Sep.py                         | 149 -------------
 dart/tests/bli_Trans.py                       | 148 -------------
 dart/tests/bli_lowLift.py                     | 195 ------------------
 dart/tests/bli_lowLift2.py                    | 170 ---------------
 dart/tests/blipython.py                       | 157 --------------
 vii/CMakeLists.txt                            |  28 +++
 vii/__init__.py                               |  21 ++
 vii/_src/CMakeLists.txt                       |  35 ++++
 vii/_src/viiw.h                               |  24 +++
 vii/_src/viiw.i                               |  53 +++++
 vii/pyVII/DAirfoil.py                         | 140 +++++++++++++
 vii/pyVII/DBoundary.py                        |  36 ++++
 vii/pyVII/DGroupConstructor.py                | 180 ++++++++++++++++
 vii/pyVII/DWake.py                            |  77 +++++++
 {dart => vii}/pyVII/vCoupler.py               |  53 ++---
 {dart => vii}/pyVII/vCoupler2.py              |   0
 {dart => vii}/pyVII/vInterpolator.py          |   0
 {dart => vii}/pyVII/vRbf.py                   |   0
 {dart => vii}/pyVII/vUtils.py                 |   8 +-
 vii/src/CMakeLists.txt                        |  27 +++
 vii/src/vii.h                                 |  50 +++++
 {dart => vii}/src/wBLRegion.cpp               |   2 +-
 {dart => vii}/src/wBLRegion.h                 |   4 +-
 {dart => vii}/src/wClosures.cpp               |   2 +-
 {dart => vii}/src/wClosures.h                 |   6 +-
 vii/src/wDartCoupler.cpp                      | 172 +++++++++++++++
 vii/src/wDartCoupler.h                        |  35 ++++
 {dart => vii}/src/wInvBnd.cpp                 |   4 +-
 {dart => vii}/src/wInvBnd.h                   |   6 +-
 {dart => vii}/src/wTimeSolver.cpp             |   2 +-
 {dart => vii}/src/wTimeSolver.h               |   4 +-
 {dart => vii}/src/wViscFluxes.cpp             |   2 +-
 {dart => vii}/src/wViscFluxes.h               |   6 +-
 {dart => vii}/src/wViscMesh.cpp               |   4 +-
 {dart => vii}/src/wViscMesh.h                 |   4 +-
 {dart => vii}/src/wViscSolver.cpp             |   4 +-
 {dart => vii}/src/wViscSolver.h               |   4 +-
 .../bli_dualMesh.py => vii/tests/bli2.py      |  82 ++++----
 {dart => vii}/tests/bli3.py                   |   0
 46 files changed, 968 insertions(+), 1275 deletions(-)
 delete mode 100644 dart/src/wInterpolator.cpp
 delete mode 100644 dart/src/wInterpolator.h
 delete mode 100644 dart/tests/bli2.py
 delete mode 100644 dart/tests/bli_Polar.py
 delete mode 100644 dart/tests/bli_Sep.py
 delete mode 100644 dart/tests/bli_Trans.py
 delete mode 100644 dart/tests/bli_lowLift.py
 delete mode 100644 dart/tests/bli_lowLift2.py
 delete mode 100755 dart/tests/blipython.py
 create mode 100644 vii/CMakeLists.txt
 create mode 100644 vii/__init__.py
 create mode 100644 vii/_src/CMakeLists.txt
 create mode 100644 vii/_src/viiw.h
 create mode 100644 vii/_src/viiw.i
 create mode 100755 vii/pyVII/DAirfoil.py
 create mode 100755 vii/pyVII/DBoundary.py
 create mode 100755 vii/pyVII/DGroupConstructor.py
 create mode 100755 vii/pyVII/DWake.py
 rename {dart => vii}/pyVII/vCoupler.py (89%)
 rename {dart => vii}/pyVII/vCoupler2.py (100%)
 rename {dart => vii}/pyVII/vInterpolator.py (100%)
 rename {dart => vii}/pyVII/vRbf.py (100%)
 rename {dart => vii}/pyVII/vUtils.py (96%)
 create mode 100644 vii/src/CMakeLists.txt
 create mode 100644 vii/src/vii.h
 rename {dart => vii}/src/wBLRegion.cpp (99%)
 rename {dart => vii}/src/wBLRegion.h (98%)
 rename {dart => vii}/src/wClosures.cpp (99%)
 rename {dart => vii}/src/wClosures.h (88%)
 create mode 100644 vii/src/wDartCoupler.cpp
 create mode 100644 vii/src/wDartCoupler.h
 rename {dart => vii}/src/wInvBnd.cpp (93%)
 rename {dart => vii}/src/wInvBnd.h (95%)
 rename {dart => vii}/src/wTimeSolver.cpp (99%)
 rename {dart => vii}/src/wTimeSolver.h (97%)
 rename {dart => vii}/src/wViscFluxes.cpp (99%)
 rename {dart => vii}/src/wViscFluxes.h (91%)
 rename {dart => vii}/src/wViscMesh.cpp (95%)
 rename {dart => vii}/src/wViscMesh.h (97%)
 rename {dart => vii}/src/wViscSolver.cpp (99%)
 rename {dart => vii}/src/wViscSolver.h (98%)
 rename dart/tests/bli_dualMesh.py => vii/tests/bli2.py (62%)
 rename {dart => vii}/tests/bli3.py (100%)

diff --git a/dart/_src/dartw.h b/dart/_src/dartw.h
index 5f7c560..54d3a31 100644
--- a/dart/_src/dartw.h
+++ b/dart/_src/dartw.h
@@ -26,11 +26,3 @@
 #include "wProblem.h"
 #include "wSolver.h"
 #include "wWake.h"
-#include "wViscSolver.h"
-#include "wBLRegion.h"
-#include "wTimeSolver.h"
-#include "wViscFluxes.h"
-#include "wClosures.h"
-#include "wInterpolator.h"
-#include "wViscMesh.h"
-#include "wInvBnd.h"
diff --git a/dart/_src/dartw.i b/dart/_src/dartw.i
index f0583cf..6847f36 100644
--- a/dart/_src/dartw.i
+++ b/dart/_src/dartw.i
@@ -58,15 +58,6 @@ threads="1"
 %shared_ptr(dart::Newton);
 %shared_ptr(dart::Adjoint);
 
-%include "wViscSolver.h"
-%include "wBLRegion.h"
-%include "wTimeSolver.h"
-%include "wViscFluxes.h"
-%include "wClosures.h"
-%include "wInterpolator.h"
-%include "wViscMesh.h"
-%include "wInvBnd.h"
-
 %include "wFluid.h"
 %include "wAssign.h"
 %include "wFreestream.h"
diff --git a/dart/src/dart.h b/dart/src/dart.h
index f242f89..ad5628c 100644
--- a/dart/src/dart.h
+++ b/dart/src/dart.h
@@ -65,16 +65,6 @@ class Newton;
 class FpeLSFunction;
 class Adjoint;
 
-// Viscous
-class ViscSolver;
-class BLRegion;
-class TimeSolver;
-class ViscSolver;
-class Closures;
-class Interpolator;
-class ViscMesh;
-class InvBnd;
-
 // General
 class F0Ps;
 class F0El;
diff --git a/dart/src/wInterpolator.cpp b/dart/src/wInterpolator.cpp
deleted file mode 100644
index afef795..0000000
--- a/dart/src/wInterpolator.cpp
+++ /dev/null
@@ -1,18 +0,0 @@
-#include "wInterpolator.h"
-#include "wClosures.h"
-#include "wViscMesh.h"
-#include "wInvBnd.h"
-#include <iostream>
-#include <vector>
-#include <cmath>
-
-using namespace dart;
-
-Interpolator::Interpolator()
-{
-} // end Interpolator
-
-Interpolator::~Interpolator()
-{
-} // end ~Interpolator
-
diff --git a/dart/src/wInterpolator.h b/dart/src/wInterpolator.h
deleted file mode 100644
index 95453bb..0000000
--- a/dart/src/wInterpolator.h
+++ /dev/null
@@ -1,22 +0,0 @@
-#ifndef WINTERPOLATOR_H
-#define WINTERPOLATOR_H
-
-#include "dart.h"
-#include "wClosures.h"
-#include "wViscMesh.h"
-#include <vector>
-#include <string>
-
-namespace dart
-{
-
-class DART_API Interpolator
-{
-    ViscMesh *vmesh;
-
-public:
-    Interpolator();
-    ~Interpolator();
-};
-} // namespace dart
-#endif //WINTERPOLATOR_H
\ No newline at end of file
diff --git a/dart/tests/bli2.py b/dart/tests/bli2.py
deleted file mode 100644
index e16ea38..0000000
--- a/dart/tests/bli2.py
+++ /dev/null
@@ -1,158 +0,0 @@
-#!/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.
-
-
-## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import dart.utils as floU
-import dart.default as floD
-#import dart.viscUtils as viscU
-
-import dart.pyVII.vCoupler as floVC
-import dart.pyVII.vUtils as viscU
-import numpy as np
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-from matplotlib import pyplot as plt
-
-import math
-
-def main():
-    print('Start')
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 1e7
-    alpha = 2*math.pi/180
-    M_inf = 0.73
-
-    CFL0 = 1
-
-    nSections = 1
-    
-    # ========================================================================================
-
-    c_ref = 1
-    dim = 2
-
-    # 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'])
-    tms['msh'].stop()
-
-    # set the problem and add medium, initial/boundary conditions
-    tms['pre'].start()
-    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-    tms['pre'].stop()
-
-    # solve viscous problem
-
-    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-    tms['solver'].start()
-    
-    isolver = floD.newton(pbl)
-    vsolver = viscU.initBL(Re, M_inf, CFL0, nSections)
-    config = {'nDim': dim, 'nSections':nSections, 'nPoints':[600, 50], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 10}
-
-    coupler = floVC.Coupler(isolver, vsolver)
-    coupler.CreateProblem(blw[0], blw[1])
-    coupler.Run()
-
-    # extract Cp
-    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
-    tboxU.write(Cp, 'Cp.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, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
-
-    # Write results to file
-    vSol = viscU.GetSolution(0, vsolver)
-    plt.plot(vSol['x/c'], vSol['Cf'])
-    plt.show()
-    viscU.WriteFile(vSol, vsolver.GetRe())
-
-    # display timers
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    print('SOLVERS statistics')
-    print(coupler.tms)
-
-    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)
-
-    blScalars, blVectors = viscU.ExtractVars(vsolver)
-    wData = viscU.Dictionary(blScalars, blVectors, xtr)
-    viscU.PlotVars(wData, plotVar)"""
-
-    
-    # check results
-    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', 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', 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()
-
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/dart/tests/bli_Polar.py b/dart/tests/bli_Polar.py
deleted file mode 100644
index 3e08f13..0000000
--- a/dart/tests/bli_Polar.py
+++ /dev/null
@@ -1,124 +0,0 @@
-#!/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.
-
-
-## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import dart.utils as floU
-import dart.default as floD
-import dart.VII.viscUtils as viscU
-
-import dart.VII.vCoupler as floVC
-import numpy as np
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-import math
-
-def main():
-    print('Start')
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 1e7
-    alpha=0
-    alphaSeq = np.arange(-5,5.5,5)
-    M_inf = 0.
-
-    meshFactor = 1
-    CFL0 = 1
-
-    # ========================================================================================
-
-    c_ref = 1
-    dim = 2
-
-    # mesh the geometry
-    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
-    tms['msh'].start()
-    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.005, 'msLe' : 0.001}
-    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
-    tms['msh'].stop()
-
-    # set the problem and add medium, initial/boundary conditions
-    tms['pre'].start()
-    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-    tms['pre'].stop()
-
-    # solve viscous problem
-
-    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-    tms['solver'].start()
-    
-    isolver = floD.newton(pbl)
-    vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
-
-    coupler = floVC.Coupler(isolver, vsolver)
-
-    coupler.CreateProblem(blw[0], blw[1])
-
-    polar = coupler.RunPolar(alphaSeq)
-    #coupler.Run()
-    tms['solver'].stop()
-
-    # display timers
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    print('SOLVERS statistics')
-    print(coupler.tms)
-
-    #viscU.PlotPolar(polar)
-    
-    # check results
-    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
-    tests = CTests()
-    tests.add(CTest('Cl1', polar[1,0], -0.430901, 5e-3))
-    tests.add(CTest('Cd1', polar[2,0], 0.005728, 1e-3, forceabs=True))
-    tests.add(CTest('Cl2', polar[1,-1], 0.430797, 5e-3))
-    tests.add(CTest('Cd2', polar[2,-1], 0.005789, 1e-3, forceabs=True))
-
-    tests.run()
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/dart/tests/bli_Sep.py b/dart/tests/bli_Sep.py
deleted file mode 100644
index cd35876..0000000
--- a/dart/tests/bli_Sep.py
+++ /dev/null
@@ -1,149 +0,0 @@
-#!/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.
-
-
-## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import dart.utils as floU
-import dart.default as floD
-import dart.VII.viscUtils as viscU
-
-import dart.VII.vCoupler as floVC
-import numpy as np
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-import math
-
-import math
-
-def main():
-    print('Start')
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 1e7
-    alpha = 15*math.pi/180
-    #alphaSeq = np.arange(-5, 5.5, 0.5)
-    M_inf = 0.2
-
-    meshFactor = 2
-    CFL0 = 1
-
-    plotVar = 'H'
-
-    # ========================================================================================
-
-    c_ref = 1
-    dim = 2
-
-    # 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.0001}
-    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
-    tms['msh'].stop()
-
-    # set the problem and add medium, initial/boundary conditions
-    tms['pre'].start()
-    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-    tms['pre'].stop()
-
-    # solve viscous problem
-
-    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-    tms['solver'].start()
-    
-    isolver = floD.newton(pbl)
-    vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
-
-    coupler = floVC.Coupler(isolver, vsolver)
-
-    coupler.CreateProblem(blw[0], blw[1])
-
-    #coupler.RunPolar(alphaSeq)
-    coupler.Run()
-    tms['solver'].stop()
-
-    # extract Cp
-    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, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
-
-    # display timers
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    
-    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)
-
-    blScalars, blVectors = viscU.ExtractVars(vsolver)
-    xtr=vsolver.Getxtr()
-    wData = viscU.Dictionary(blScalars, blVectors, xtr)
-    viscU.PlotVars(wData, plotVar)"""
-
-    
-    # check results
-    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
-    tests = CTests()
-    if Re == 1e7 and M_inf == 0.2 and alpha == 15*math.pi/180:
-        tests.add(CTest('Cl', isolver.Cl, 1.61416, 5e-2))
-        tests.add(CTest('Cd', vsolver.Cdt, 0.01452, 1e-3, forceabs=True))
-        tests.add(CTest('Cdp', vsolver.Cdp, 0.0105, 1e-3, forceabs=True))
-        tests.add(CTest('xtr_top', xtr[0], 0.01, 0.03, forceabs=True))
-        tests.add(CTest('xtr_bot', xtr[1], 1.00, 0.03, forceabs=True))
-    else:
-        raise Exception('Test not defined for this flow')
-
-    tests.run()
-
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/dart/tests/bli_Trans.py b/dart/tests/bli_Trans.py
deleted file mode 100644
index 9e27a0a..0000000
--- a/dart/tests/bli_Trans.py
+++ /dev/null
@@ -1,148 +0,0 @@
-#!/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.
-
-
-## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import dart.utils as floU
-import dart.default as floD
-import dart.VII.viscUtils as viscU
-
-import dart.VII.vCoupler as floVC
-import numpy as np
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-import math
-
-def main():
-    print('Start')
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 1e7
-    alpha = 2.*math.pi/180
-    #alphaSeq = np.arange(-5, 5.5, 0.5)
-    M_inf = 0.715
-
-    meshFactor = 2
-    CFL0 = 1
-
-    plotVar = 'H'
-
-    # ========================================================================================
-
-    c_ref = 1
-    dim = 2
-
-    # 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'])
-    tms['msh'].stop()
-
-    # set the problem and add medium, initial/boundary conditions
-    tms['pre'].start()
-    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-    tms['pre'].stop()
-
-    # solve viscous problem
-
-    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-    tms['solver'].start()
-    
-    isolver = floD.newton(pbl)
-    vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
-
-    coupler = floVC.Coupler(isolver, vsolver)
-
-    coupler.CreateProblem(blw[0], blw[1])
-
-    #coupler.RunPolar(alphaSeq)
-    coupler.resetInv = True
-    coupler.Run()
-    tms['solver'].stop()
-
-    # extract Cp
-    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, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
-
-    # display timers
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    
-    
-    
-    """# 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)
-
-    blScalars, blVectors = viscU.ExtractVars(vsolver)
-    xtr=vsolver.Getxtr()
-    wData = viscU.Dictionary(blScalars, blVectors, xtr)
-    viscU.PlotVars(wData, plotVar)"""
-
-    
-    # check results
-    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
-    tests = CTests()
-    if Re == 1e7 and M_inf == 0.2 and alpha == 5*math.pi/180:
-        tests.add(CTest('Cl', isolver.Cl, 0.56, 5e-2)) # XFOIL 0.58
-        tests.add(CTest('Cd', vsolver.Cdt, 0.0063, 1e-3, forceabs=True)) # XFOIL 0.00617
-        tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True)) # XFOIL 0.00176
-        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.059, 0.03, forceabs=True)) # XFOIL 0.0510
-        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.738, 0.03, forceabs=True)) # XFOIL 0.7473
-    else:
-        raise Exception('Test not defined for this flow')
-
-    tests.run()
-
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
deleted file mode 100644
index 10c1ae6..0000000
--- a/dart/tests/bli_lowLift.py
+++ /dev/null
@@ -1,195 +0,0 @@
-#!/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.
-
-
-## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import dart.utils as floU
-import dart.default as floD
-#import dart.viscUtils as viscU
-
-import dart.pyVII.vCoupler as floVC
-import dart.pyVII.vUtils as viscU
-import numpy as np
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-from matplotlib import pyplot as plt
-
-import math
-
-def main():
-    print('Start')
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 1e7
-    alpha = 5*math.pi/180
-    M_inf = 0.
-
-    CFL0 = 1
-
-    nSections = 1
-    
-    # ========================================================================================
-
-    c_ref = 1
-    dim = 2
-
-    nRel = 1
-
-    nomDeVariableDePorc_x = [None for _ in range(nRel)]
-    nomDeVariableDePorc_var = [None for _ in range(nRel)]
-    nomDeVariableDePorc_Uin = [None for _ in range(nRel)]
-
-    for i in range(nRel):
-        # 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.001}
-        msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
-        tms['msh'].stop()
-
-        # set the problem and add medium, initial/boundary conditions
-        tms['pre'].start()
-        pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-        tms['pre'].stop()
-
-        # solve viscous problem
-
-        print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-        tms['solver'].start()
-        
-        isolver = floD.newton(pbl)
-        vsolver = viscU.initBL(Re, M_inf, CFL0, nSections)
-        config = {'nDim': dim, 'nSections':nSections, 'nPoints':[600, 50], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 10}
-
-        coupler = floVC.Coupler(isolver, vsolver)
-        coupler.CreateProblem(blw[0], blw[1])
-        #fig, axs = plt.subplots(2)
-        aeroCoeff = coupler.Run()
-        #isolver.reset()
-
-        nomDeVariableDePorc_x[i] = coupler.sln[15].copy()
-        nomDeVariableDePorc_var[i] = coupler.sln[1].copy()
-        nomDeVariableDePorc_Uin[i] = coupler.fMeshDict['vx'].copy()
-
-        del isolver
-        del vsolver
-        del coupler
-        del pbl
-        del bnd
-        del blw
-        del msh
-        del gmshWriter
-
-    print(len(nomDeVariableDePorc_x), len(nomDeVariableDePorc_var))
-    for i in range(len(nomDeVariableDePorc_Uin)):
-        plt.plot(nomDeVariableDePorc_Uin[i] - nomDeVariableDePorc_Uin[0], 'x-')
-    plt.show()
-    for i in range(len(nomDeVariableDePorc_var)):
-        plt.plot(nomDeVariableDePorc_x[i], nomDeVariableDePorc_var[i] - nomDeVariableDePorc_var[0], 'x-')
-    plt.show()
-
-    for i in range(len(nomDeVariableDePorc_var)):
-        plt.plot(nomDeVariableDePorc_x[i], nomDeVariableDePorc_var[i], '-')
-    plt.show()
-
-    """sln = np.asarray(vsolver.GetSolution(0))
-
-    print(sln[8][sln[8]>0.1])
-    plt.plot(sln[16], sln[1])
-    plt.xlim([0, 2])
-    plt.show()"""
-    #coupler.RunPolar(alphaSeq)
-    tms['solver'].stop()
-
-    # extract Cp
-    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, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
-
-    # display timers
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    print('SOLVERS statistics')
-    print(coupler.tms)
-
-    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)
-
-    blScalars, blVectors = viscU.ExtractVars(vsolver)
-    wData = viscU.Dictionary(blScalars, blVectors, xtr)
-    viscU.PlotVars(wData, plotVar)"""
-
-    
-    # check results
-    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', 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', 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()
-
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/dart/tests/bli_lowLift2.py b/dart/tests/bli_lowLift2.py
deleted file mode 100644
index 05ab0ab..0000000
--- a/dart/tests/bli_lowLift2.py
+++ /dev/null
@@ -1,170 +0,0 @@
-#!/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.
-
-
-## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-from email.errors import NonPrintableDefect
-import dart.utils as floU
-import dart.default as floD
-#import dart.viscUtils as viscU
-
-import dart.pyVII.vCoupler2 as viscC # Rbf coupler
-import dart.pyVII.vInterpolator as Interpol # Rbf interpolator
-import dart.pyVII.vUtils as viscU
-import numpy as np
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-from matplotlib import pyplot as plt
-
-import math
-
-def main():
-    print('Start')
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 1e7
-    alpha = 5*math.pi/180
-    M_inf = 0.
-
-    CFL0 = 1
-
-    nSections = 1
-    verbose = 0
-    
-    # ========================================================================================
-
-    c_ref = 1
-    dim = 2
-
-    nPv = [200]
-    aeroCoeff = [np.zeros(2) for i in range(len(nPv))]
-
-
-    for i in range(len(nPv)):
-
-        # 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.001}
-        msh, sortElm, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
-        tms['msh'].stop()
-
-        # set the problem and add medium, initial/boundary conditions
-        tms['pre'].start()
-        pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-        tms['pre'].stop()
-
-        # solve viscous problem
-
-        print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-        tms['solver'].start()
-        
-        isolver = floD.newton(pbl)
-        if verbose == 0:
-            isolver.printOn = False
-        vsolver = viscU.initBL(Re, M_inf, CFL0, nSections, verbose)
-        config = {'nDim': dim, 'nSections':nSections, 'nPoints':[nPv[i], 25], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 15}
-
-        iSolverAPI = Interpol.Interpolator(isolver, sortElm, blw, config)
-
-        coupler = viscC.Coupler(iSolverAPI, vsolver)
-        aeroCoeff[i] = coupler.Run()
-
-        """del isolver
-        del vsolver
-        del iSolverAPI
-        del pbl
-        del bnd
-        del blw
-        del msh
-        del gmshWriter"""
-    tms['solver'].stop()
-
-    # extract Cp
-    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, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
-
-    # display timers
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    print('SOLVERS statistics')
-    print(coupler.tms)
-
-    xtr=[vsolver.Getxtr(0,0),vsolver.Getxtr(0,1)] 
-    
-    # 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)
-
-    
-    # check results
-    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', 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', 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()
-
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/dart/tests/blipython.py b/dart/tests/blipython.py
deleted file mode 100755
index 8f40e54..0000000
--- a/dart/tests/blipython.py
+++ /dev/null
@@ -1,157 +0,0 @@
-#!/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.
-
-
-## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import dart.utils as floU
-import dart.default as floD
-
-import dart.viscous.Master.DCoupler as floVC
-import dart.viscous.Drivers.DDriver as floVSMaster
-
-import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
-import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
-
-import dart.viscous.Physics.DValidation as validation
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-import math
-
-def main():
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 1e7
-    alpha = 12.*math.pi/180
-    M_inf = 0.
-
-    #user_xtr=[0.41,0.41]
-    #user_xtr=[0,None]
-    user_xtr=[None,None]
-    user_Ti=None
-
-    mapMesh = 0
-    nFactor = 2
-
-    # Time solver parameters
-    Time_Domain = False # True for unsteady solver, False for steady solver
-    CFL0 = 1
-
-    # ========================================================================================
-
-    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
-    c_ref = 1
-    dim = 2
-    tol = 1e-4 # tolerance for the VII (usually one drag count)
-    case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
-
-    # 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}
-    #pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.002, 'msLe' : 0.0001}
-    msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
-    tms['msh'].stop()
-
-    # set the problem and add medium, initial/boundary conditions
-    tms['pre'].start()
-    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', vsc = True, dbc=True)
-    tms['pre'].stop()
-
-    # solve viscous problem
-
-    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-    tms['solver'].start()
-    isolver = floD.newton(pbl)
-
-    # Choose between steady and unsteady solver
-    if Time_Domain is True: # Run unsteady solver
-        vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
-        coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
-
-    elif Time_Domain is False: # Run steady solver 
-        vsolver=floVS_Std.Solver(Re)
-        coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
-    coupler.run()
-    tms['solver'].stop()
-
-    # extract Cp
-    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, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
-
-    # display timers
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    print('SOLVERS statistics')
-    print(coupler.tms)
-    
-    # check results
-    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', isolver.Cl, 0.2197, 5e-2)) # XFOIL 0.58
-        tests.add(CTest('Cd', vsolver.Cd, 0.00529, 1e-3, forceabs=True)) # XFOIL 0.00617
-        tests.add(CTest('Cdp', vsolver.Cdp, 0.0008, 1e-3, forceabs=True)) # XFOIL 0.00176
-        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.20395, 0.03, forceabs=True)) # XFOIL 0.0510
-        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.49996, 0.03, forceabs=True)) # XFOIL 0.7473
-    else:
-        raise Exception('Test not defined for this flow')
-
-    tests.run()
-
-    # 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.Cd, isolver.Cm, 4), True)
-
-    val=validation.Validation(case)
-    val.main(isolver.Cl,vsolver.wData)
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/vii/CMakeLists.txt b/vii/CMakeLists.txt
new file mode 100644
index 0000000..8de8317
--- /dev/null
+++ b/vii/CMakeLists.txt
@@ -0,0 +1,28 @@
+# 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.
+
+# Add source dir
+ADD_SUBDIRECTORY( src )
+ADD_SUBDIRECTORY( _src )
+
+# Add test dir
+MACRO_AddTest(${CMAKE_CURRENT_SOURCE_DIR}/tests)
+
+# Add to install
+INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py
+              ${CMAKE_CURRENT_LIST_DIR}/utils.py
+        DESTINATION vii)
+INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/api
+                  ${CMAKE_CURRENT_LIST_DIR}/viscous
+        DESTINATION vii)
diff --git a/vii/__init__.py b/vii/__init__.py
new file mode 100644
index 0000000..43be8da
--- /dev/null
+++ b/vii/__init__.py
@@ -0,0 +1,21 @@
+# -*- coding: utf-8 -*-
+# dart MODULE initialization file
+
+# 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.
+
+
+import fwk
+import tbox
+from viiw import *
diff --git a/vii/_src/CMakeLists.txt b/vii/_src/CMakeLists.txt
new file mode 100644
index 0000000..18dc158
--- /dev/null
+++ b/vii/_src/CMakeLists.txt
@@ -0,0 +1,35 @@
+# 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.
+
+# CMake input file of the SWIG wrapper around "viiw.so"
+
+INCLUDE(${SWIG_USE_FILE})
+
+FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
+FILE(GLOB ISRCS *.i)
+
+SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
+SET(CMAKE_SWIG_FLAGS "-interface" "_viiw") # avoids "import _module_d" with MSVC/Debug
+SWIG_ADD_LIBRARY(viiw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
+SET_PROPERTY(TARGET viiw PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON)
+MACRO_DebugPostfix(viiw)
+
+TARGET_INCLUDE_DIRECTORIES(viiw PRIVATE ${PROJECT_SOURCE_DIR}/vii/_src
+                                         ${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
+                                         ${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
+                                         ${PYTHON_INCLUDE_PATH})
+TARGET_LINK_LIBRARIES(viiw PRIVATE vii tbox fwk ${PYTHON_LIBRARIES})
+
+INSTALL(FILES ${EXECUTABLE_OUTPUT_PATH}/\${BUILD_TYPE}/viiw.py DESTINATION ${CMAKE_INSTALL_PREFIX})
+INSTALL(TARGETS viiw DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/vii/_src/viiw.h b/vii/_src/viiw.h
new file mode 100644
index 0000000..8060687
--- /dev/null
+++ b/vii/_src/viiw.h
@@ -0,0 +1,24 @@
+/*
+ * 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.
+ */
+
+#include "wViscSolver.h"
+#include "wBLRegion.h"
+#include "wTimeSolver.h"
+#include "wViscFluxes.h"
+#include "wClosures.h"
+#include "wViscMesh.h"
+#include "wInvBnd.h"
+#include "wCoupler.h"
\ No newline at end of file
diff --git a/vii/_src/viiw.i b/vii/_src/viiw.i
new file mode 100644
index 0000000..2d788b1
--- /dev/null
+++ b/vii/_src/viiw.i
@@ -0,0 +1,53 @@
+/*
+ * 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.
+ */
+
+// SWIG input file of the 'dart' module
+
+%feature("autodoc","1");
+
+%module(docstring=
+"'viiw' module: Viscous inviscid interaction for dart 2021-2022
+(c) ULiege - A&M",
+directors="0",
+threads="1"
+) viiw
+%{
+
+#include "vii.h"
+
+#include "fwkw.h"
+#include "tboxw.h"
+#include "viiw.h"
+
+%}
+
+
+%include "fwkw.swg"
+
+// ----------- MODULES UTILISES ------------
+%import "tboxw.i"
+
+// ----------- FLOW CLASSES ----------------
+%include "vii.h"
+
+%include "wViscSolver.h"
+%include "wBLRegion.h"
+%include "wTimeSolver.h"
+%include "wViscFluxes.h"
+%include "wClosures.h"
+%include "wViscMesh.h"
+%include "wInvBnd.h"
+%include "wCoupler.h"
\ No newline at end of file
diff --git a/vii/pyVII/DAirfoil.py b/vii/pyVII/DAirfoil.py
new file mode 100755
index 0000000..e8986c8
--- /dev/null
+++ b/vii/pyVII/DAirfoil.py
@@ -0,0 +1,140 @@
+#!/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.
+
+
+## Airfoil around which the boundary layer is computed
+# Amaury Bilocq
+
+from vii.pyVII.DBoundary import Boundary
+
+import numpy as np
+
+class Airfoil(Boundary):
+    def __init__(self, _boundary):
+        Boundary.__init__(self, _boundary)
+        self.name = 'airfoil' # type of boundary
+        self.T0 = 0 # initial condition for the momentum thickness
+        self.H0 = 0
+        self.n0 = 0
+        self.Ct0 = 0
+        self.TEnd = [0,0]
+        self.HEnd = [0,0]
+        self.CtEnd = [0,0]
+        self.nEnd = [0,0]
+    
+    def initialConditions(self, Re, dv):
+        if dv > 0:
+            self.T0 = np.sqrt(0.075/(Re*dv))
+        else:
+            self.T0 = 1e-6
+        self.H0 = 2.23 # One parameter family Falkner Skan
+        self.n0 = 0
+        self.Ct0 = 0
+        return self.T0, self.H0, self.n0, self.Ct0
+
+
+    def connectList(self, couplIter):
+        ''' Sort the value read by the viscous solver/ Create list of connectivity
+        '''
+        N1 = np.zeros(self.nN, dtype=int) # Node number
+        connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
+        connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
+        data = np.zeros((self.boundary.nodes.size(), 10))
+        i = 0
+        for n in self.boundary.nodes:
+            data[i,0] = n.no
+            data[i,1] = n.pos[0]
+            data[i,2] = n.pos[1]
+            data[i,3] = n.pos[2]
+            data[i,4] = self.v[i,0]
+            data[i,5] = self.v[i,1]
+            data[i,6] = self.Me[i]
+            data[i,7] = self.rhoe[i]
+            data[i,8] = self.deltaStar[i]
+            data[i,9] = self.xx[i]
+            i += 1
+        # Table containing the element and its nodes
+        eData = np.zeros((self.nE,3), dtype=int)
+        for i in range(0, self.nE):
+            eData[i,0] = self.boundary.tag.elems[i].no
+            eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
+            eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
+        # Find the stagnation point
+        if couplIter == 0:
+          self.idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
+          idxStag = self.idxStag
+        else:
+          idxStag = self.idxStag
+        #idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
+        globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
+        # Find TE nodes (assuming suction side element is above pressure side element in the y direction)
+        idxTE = np.where(data[:,1] == np.max(data[:,1]))[0] # TE nodes
+        idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], eData[:,2] == data[idxTE[0],0]))[0] # TE element 1
+        idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], eData[:,2] == data[idxTE[1],0]))[0] # TE element 2
+        ye0 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe0[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe0[0],2])[0],2]) # y coordinates of TE element 1 CG
+        ye1 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe1[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe1[0],2])[0],2]) # y coordinates of TE element 2 CG
+        if ye0 - ye1 > 0: # element 1 containing TE node 1 is above element 2
+            upperGlobTE = int(data[idxTE[0],0])
+            lowerGlobTE = int(data[idxTE[1],0])
+        else:
+             upperGlobTE = int(data[idxTE[1],0])
+             lowerGlobTE = int(data[idxTE[0],0])
+        # Connectivity
+        connectListElems[0] = np.where(eData[:,1] == globStag)[0]
+        N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes -> globStag
+        connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
+        i = 1
+        upperTE = 0
+        lowerTE = 0
+        # Sort the suction part
+        while upperTE == 0:
+            N1[i] = eData[connectListElems[i-1],2] # Second node of the element
+            connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes
+            connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
+            if eData[connectListElems[i],2] == upperGlobTE:
+                upperTE = 1
+            i += 1
+        # Sort the pressure side
+        connectListElems[i] = np.where(eData[:,2] == globStag)[0]
+        connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
+        N1[i] = eData[connectListElems[i],2]
+        while lowerTE == 0:
+            N1[i+1]  = eData[connectListElems[i],1] # First node of the element
+            connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes
+            connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
+            if eData[connectListElems[i+1],1] == lowerGlobTE:
+                lowerTE = 1
+            i += 1
+        connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
+        data[:,0] = data[connectListNodes,0]
+        data[:,1] = data[connectListNodes,1]
+        data[:,2] = data[connectListNodes,2]
+        data[:,3] = data[connectListNodes,3]
+        data[:,4] = data[connectListNodes,4]
+        data[:,5] = data[connectListNodes,5]
+        data[:,6] = data[connectListNodes,6]
+        data[:,7] = data[connectListNodes,7]
+        data[:,8] = data[connectListNodes,8]
+        data[:,9] = data[connectListNodes,9]
+        # Separated upper and lower part
+        data = np.delete(data,0,1)
+        uData = data[0:np.argmax(data[:,0])+1]
+        lData = data[np.argmax(data[:,0])+1:None]
+        lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
+        print(np.shape(connectListNodes))
+        quit()
+        return connectListNodes, connectListElems,[uData, lData]
diff --git a/vii/pyVII/DBoundary.py b/vii/pyVII/DBoundary.py
new file mode 100755
index 0000000..f3ff33a
--- /dev/null
+++ b/vii/pyVII/DBoundary.py
@@ -0,0 +1,36 @@
+#!/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.
+
+
+## Base class representing a physical boundary
+# Amaury Bilocq
+
+import numpy as np
+
+class Boundary:
+    def __init__(self, _boundary):
+        self.boundary = _boundary # boundary of interest
+        self.nN = len(self.boundary.nodes) # number of nodes
+        self.nE = len(self.boundary.tag.elems) # number of elements
+        self.v = np.zeros((self.nN, 3)) # velocity at edge of BL
+        self.u = np.zeros(self.nE) # blowing velocity
+        self.Me = np.zeros(self.nN) # Mach number at the edge of the boundary layer
+        self.rhoe = np.zeros(self.nN) # density at the edge of the boundary layer
+        self.connectListnodes = np.zeros(self.nN) # connectivity for nodes
+        self.connectListElems = np.zeros(self.nE) # connectivity for elements
+        self.deltaStar = np.zeros(self.nN) # displacement thickness
+        self.xx = np.zeros(self.nN) # coordinates in the reference of the physical surface. Used to store the values for the interaction law
diff --git a/vii/pyVII/DGroupConstructor.py b/vii/pyVII/DGroupConstructor.py
new file mode 100755
index 0000000..036d604
--- /dev/null
+++ b/vii/pyVII/DGroupConstructor.py
@@ -0,0 +1,180 @@
+from matplotlib import pyplot as plt
+import numpy as np
+
+class GroupConstructor():
+    """ ---------- Group Constructor ----------
+        
+        Builds single vectors elements out of the 3 groups given in input.
+    """
+    def __init__(self) -> None:
+        pass
+
+    def mergeVectors(self,dataIn, mapMesh, nFactor):
+        """Merges 3 groups upper/lower sides and wake.
+        """
+        
+        for k in range(len(dataIn)):
+
+            if len(dataIn[k]) == 2: # Airfoil case.
+            
+                xx_up, dv_up, vtExt_up, _, alpha_up = self.__getBLcoordinates(dataIn[k][0])
+                xx_lw, dv_lw, vtExt_lw, _, alpha_lw = self.__getBLcoordinates(dataIn[k][1])
+
+            else: # Wake case
+
+                xx_wk, dv_wk, vtExt_wk, _, alpha_wk = self.__getBLcoordinates(dataIn[k])
+
+        if mapMesh:
+            # Save initial mesh.
+            cMesh        = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+            cMeshbNodes  = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
+            cxx          = np.concatenate((xx_up, xx_lw[1:], xx_wk))
+            cvtExt       = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
+            cxxExt       = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][1:,8],dataIn[1][:,8]))
+
+            # Create fine mesh.
+            fMeshUp      = self.createFineMesh(dataIn[0][0][:,0], nFactor)
+            fMeshLw      = self.createFineMesh(dataIn[0][1][:,0], nFactor)
+            fMeshWk      = self.createFineMesh(dataIn[1][:,0]   , nFactor)
+            fMesh        = np.concatenate((fMeshUp, fMeshLw[1:], fMeshWk))
+            fMeshbNodes  = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw[1:])]
+
+            fxx_up       = self.interpolateToFineMesh(xx_up, fMeshUp, nFactor)
+            fxx_lw       = self.interpolateToFineMesh(xx_lw, fMeshLw, nFactor)
+            fxx_wk       = self.interpolateToFineMesh(xx_wk, fMeshWk, nFactor)
+            fxx          = np.concatenate((fxx_up, fxx_lw[1:], fxx_wk))
+
+            fvtExt_up    = self.interpolateToFineMesh(vtExt_up, fMeshUp, nFactor)
+            fvtExt_lw    = self.interpolateToFineMesh(vtExt_lw, fMeshLw, nFactor)
+            fvtExt_wk    = self.interpolateToFineMesh(vtExt_wk, fMeshWk, nFactor)
+            fvtExt       = np.concatenate((fvtExt_up, fvtExt_lw[1:], fvtExt_wk))
+
+            fMe_up       = self.interpolateToFineMesh(dataIn[0][0][:,5], fMeshUp, nFactor)
+            fMe_lw       = self.interpolateToFineMesh(dataIn[0][1][:,5], fMeshLw, nFactor)
+            fMe_wk       = self.interpolateToFineMesh(dataIn[1][:,5]   , fMeshWk, nFactor)
+            fMe          = np.concatenate((fMe_up, fMe_lw[1:], fMe_wk))
+
+            frho_up      = self.interpolateToFineMesh(dataIn[0][0][:,6], fMeshUp, nFactor)
+            frho_lw      = self.interpolateToFineMesh(dataIn[0][1][:,6], fMeshLw, nFactor)
+            frho_wk      = self.interpolateToFineMesh(dataIn[1][:,6]   , fMeshWk, nFactor)
+            frho         = np.concatenate((frho_up, frho_lw[1:], frho_wk))
+
+            fdStarExt_up = self.interpolateToFineMesh(dataIn[0][0][:,7], fMeshUp, nFactor)
+            fdStarExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,7], fMeshLw, nFactor)
+            fdStarExt_wk = self.interpolateToFineMesh(dataIn[1][:,7]   , fMeshWk, nFactor)
+            
+            fxxExt_up    = self.interpolateToFineMesh(dataIn[0][0][:,8], fMeshUp, nFactor)
+            fxxExt_lw    = self.interpolateToFineMesh(dataIn[0][1][:,8], fMeshLw, nFactor)
+            fxxExt_wk    = self.interpolateToFineMesh(dataIn[1][:,8]   , fMeshWk, nFactor)
+            
+            fdStarext    = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk))
+            fxxext       = np.concatenate((fxxExt_up, fxxExt_lw[1:], fxxExt_wk))
+
+            fdv          = np.concatenate((dv_up, dv_lw[1:], dv_wk))
+
+            fMeshDict    = {'globCoord': fMesh, 'bNodes': fMeshbNodes, 'locCoord': fxx,
+                            'vtExt': fvtExt, 'Me': fMe, 'rho': frho, 'deltaStarExt': fdStarext, 'xxExt': fxxext, 'dv': fdv}
+
+            cMe          = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][1:,5], dataIn[1][:,5]))
+            cRho         = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][1:,6], dataIn[1][:,6]))
+            cdStarExt    = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][1:,7], dataIn[1][:,7]))
+            cMeshDict    = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx,
+                            'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv}
+
+            dataUpper = np.zeros((len(fMeshUp), 0))
+            dataLower = np.zeros((len(fMeshLw), 0))
+            dataWake  = np.zeros((len(fMeshWk), 0))
+            for iParam in range(len(dataIn[0][0][0,:])):
+                interpParamUp = self.interpolateToFineMesh(dataIn[0][0][:,iParam], fMeshUp, nFactor)
+                dataUpper     = np.column_stack((dataUpper, interpParamUp))
+                interpParamLw = self.interpolateToFineMesh(dataIn[0][1][:,iParam], fMeshLw, nFactor)
+                dataLower     = np.column_stack((dataLower, interpParamLw))
+                interpParamWk = self.interpolateToFineMesh(dataIn[1][:,iParam]   , fMeshWk, nFactor)
+                dataWake      = np.column_stack((dataWake, interpParamWk))
+            
+            dataLower = np.delete(dataLower, 0, axis = 0) # Remove stagnation point doublon. 
+
+        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))
+
+            Me           = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][:,5], dataIn[1][:,5]))
+            rho          = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][:,6], dataIn[1][:,6]))
+            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, '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()
+
+            dataUpper    = dataIn[0][0]
+            dataLower    = dataIn[0][1][1:, :]
+            dataWake     = dataIn[1]
+        
+        data = np.concatenate((dataUpper, dataLower, dataWake), axis = 0)
+
+        return fMeshDict, cMeshDict, data
+
+    def __getBLcoordinates(self,data):
+        """Transform the reference coordinates into airfoil coordinates
+        """
+        nN    = len(data[:,0])
+        nE    = nN-1            # Nbr of element along the surface.
+        vt    = np.zeros(nN)    # Outer flow velocity. 
+        xx    = np.zeros(nN)    # Position along the surface of the airfoil.
+        dx    = np.zeros(nE)    # Distance along two nodes.
+        dv    = np.zeros(nE)    # Speed difference according to element.
+        alpha = np.zeros(nE)    # Angle of the element for the rotation matrix.
+
+        for i in range(0, nE):
+
+            alpha[i] = np.arctan2((data[i+1, 1] - data[i, 1]), (data[i+1, 0]-data[i, 0]))
+            dx[i]    = np.sqrt((data[i+1, 0] - data[i,0])**2 + (data[i+1, 1]-data[i, 1])**2)
+            xx[i+1]  = dx[i] + xx[i]
+            vt[i]    = (data[i, 3] * np.cos(alpha[i]) + data[i ,4] * np.sin(alpha[i]))          # Tangential speed at node 1 element i
+            vt[i+1]  = (data[i+1, 3] * np.cos(alpha[i]) + data[i+1, 4] * np.sin(alpha[i]))      # Tangential speed at node 2 element i
+            dv[i]    = (vt[i+1] - vt[i]) / dx[i]                                                # Velocity gradient according to element i.
+                                                                                                #             central scheme with half point
+        
+        return xx, dv, vt, nN, alpha
+
+    def interpolateToFineMesh(self, cVector, fMesh, nFactor):
+        
+        fVector = np.zeros(len(fMesh)-1)
+        for cPoint in range(len(cVector) - 1):
+            dv = cVector[cPoint + 1] - cVector[cPoint]
+
+            for fPoint in range(nFactor):
+                fVector[nFactor * cPoint + fPoint] = cVector[cPoint] + fPoint * dv / nFactor
+        
+        fVector = np.append(fVector, cVector[-1])
+
+        return fVector
+
+    def createFineMesh(self, cMesh, nFactor):
+
+        fMesh = []
+        for iPoint in range(len(cMesh) - 1):
+
+            dx = cMesh[iPoint + 1] - cMesh[iPoint]
+
+            for iRefinedPoint in range(nFactor):
+
+                fMesh = np.append(fMesh, cMesh[iPoint] + iRefinedPoint * dx / nFactor)
+
+        fMesh = np.append(fMesh, cMesh[-1])
+        return fMesh
diff --git a/vii/pyVII/DWake.py b/vii/pyVII/DWake.py
new file mode 100755
index 0000000..d9f0540
--- /dev/null
+++ b/vii/pyVII/DWake.py
@@ -0,0 +1,77 @@
+#!/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.
+
+
+## Wake behind airfoil (around which the boundary layer is computed)
+# Amaury Bilocq
+
+from vii.pyVII.DBoundary import Boundary
+
+import numpy as np
+
+class Wake(Boundary):
+    def __init__(self, _boundary):
+        Boundary.__init__(self, _boundary)
+        self.name = 'wake' # type of boundary
+        self.T0 = 0 # initial condition for the momentum thickness
+        self.H0 = 0
+        self.n0 = 9 # wake is always turbulent
+        self.Ct0 = 0
+
+    def connectList(self, couplIter):
+        ''' Sort the value read by the viscous solver/ Create list of connectivity
+        '''
+        N1 = np.zeros(self.nN, dtype=int) # Node number
+        connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
+        connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
+        data = np.zeros((self.boundary.nodes.size(), 10))
+        i = 0
+        for n in self.boundary.nodes:
+            data[i,0] = n.no
+            data[i,1] = n.pos[0]
+            data[i,2] = n.pos[1]
+            data[i,3] = n.pos[2]
+            data[i,4] = self.v[i,0]
+            data[i,5] = self.v[i,1]
+            data[i,6] = self.Me[i]
+            data[i,7] = self.rhoe[i]
+            data[i,8] = self.deltaStar[i]
+            data[i,9] = self.xx[i]
+            i += 1
+        # Table containing the element and its nodes
+        eData = np.zeros((self.nE,3), dtype=int)
+        for i in range(0, self.nE):
+            eData[i,0] = self.boundary.tag.elems[i].no
+            eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
+            eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
+        connectListNodes = data[:,1].argsort()
+        connectListElems[0] = np.where(eData[:,1] == min(data[:,1]))[0]
+        for i in range(1, len(eData[:,0])):
+            connectListElems[i] = np.where(eData[:,1] == eData[connectListElems[i-1],2])[0]
+        data[:,0] = data[connectListNodes,0]
+        data[:,1] = data[connectListNodes,1]
+        data[:,2] = data[connectListNodes,2]
+        data[:,3] = data[connectListNodes,3]
+        data[:,4] = data[connectListNodes,4]
+        data[:,5] = data[connectListNodes,5]
+        data[:,6] = data[connectListNodes,6]
+        data[:,7] = data[connectListNodes,7]
+        data[:,8] = data[connectListNodes,8]
+        data[:,9] = data[connectListNodes,9]
+        # Separated upper and lower part
+        data = np.delete(data,0,1)
+        return connectListNodes, connectListElems, data
diff --git a/dart/pyVII/vCoupler.py b/vii/pyVII/vCoupler.py
similarity index 89%
rename from dart/pyVII/vCoupler.py
rename to vii/pyVII/vCoupler.py
index f74cd4c..78d78c9 100644
--- a/dart/pyVII/vCoupler.py
+++ b/vii/pyVII/vCoupler.py
@@ -22,9 +22,9 @@
 import dart.utils as floU
 import tbox.utils as tboxU
 
-from dart.viscous.Master.DAirfoil import Airfoil
-from dart.viscous.Master.DWake import Wake
-from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
+from vii.pyVII.DAirfoil import Airfoil
+from vii.pyVII.DWake import Wake
+from vii.pyVII.DGroupConstructor import GroupConstructor
 
 from fwk.coloring import ccolors
 import math
@@ -35,15 +35,13 @@ import numpy as np
 from matplotlib import pyplot as plt
 
 class Coupler:
-  def __init__(self, iSolver, vSolver) -> None:
+  def __init__(self, iSolver, vSolver, _iterPrint=1) -> None:
     self.iSolver=iSolver
     self.vSolver=vSolver
 
-    self.vSolver.verbose = self.iSolver.verbose
-    #self.vSolver.verbose = 0
-
     self.maxCouplIter = 150
     self.couplTol = 1e-4
+    self.iterPrint = _iterPrint
 
     self.tms=fwk.Timers()
     self.resetInv = False
@@ -61,18 +59,17 @@ class Coupler:
 
     # Initilize meshes in c++ objects
 
-    #self.iSolver.printOn = False
-    #self.vSolver.printOn = False
-
-    dragIter = np.zeros((0,2))
+    aeroCoeffs = np.zeros((0,2))
 
     self.nElmUpperPrev = 0
     couplIter = 0
     cdPrev = 0.0
-    """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('- 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.GetRe()/1e6))
+    print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
     print(' --------------------------------------------------------')
+
     while couplIter < self.maxCouplIter:
 
       # Run inviscid solver
@@ -96,10 +93,6 @@ class Coupler:
       self.tms['viscous'].start()
       exitCode = self.vSolver.Run(couplIter)
       self.tms['viscous'].stop()
-      # print('Viscous ops terminated with exit code ', exitCode)
-
-      #plt.xlim([0, 2])
-      #plt.show()
 
       self.tms['processing'].start()
       self.__ImposeBlwVel(couplIter)
@@ -107,35 +100,29 @@ class Coupler:
 
       error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
 
-      dragIter = np.vstack((dragIter, [self.iSolver.Cl, self.vSolver.Cdt]))
+      aeroCoeffs = np.vstack((aeroCoeffs, [self.iSolver.Cl, self.vSolver.Cdt]))
       
       if error  <= self.couplTol:
         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 dragIter
+        self.iSolver.Cl, self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)),ccolors.ANSI_RESET)
+        return aeroCoeffs
       
       cdPrev = self.vSolver.Cdt
 
-      if couplIter%1==0:
-        xtr = [1, 1]
-        xtr[0]=self.vSolver.Getxtr(0,0)
-        xtr[1]=self.vSolver.Getxtr(0,1)
+      if couplIter%self.iterPrint==0:
         print(' {:>4.0f}| {:>10.4f} {:>10.4f} | {:>10.4f} {:>8.4f} {:>8.3f}'.format(couplIter,
-        self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)))
+        self.iSolver.Cl, self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)))
+
+        if self.iSolver.verbose != 0 or self.vSolver.verbose != 0:
+          print('')
 
       couplIter += 1
       self.tms['processing'].stop()
-
-      self.sln = np.asarray(self.vSolver.GetSolution(0))
-      """if couplIter==1:
       
-        #plt.plot(sln[16], sln[1], 'x-')
-        #plt.plot(self.fMeshDict['globCoord'][:self.fMeshDict['bNodes'][1]], self.fMeshDict['vx'][:self.fMeshDict['bNodes'][1]], 'x-')
-        return dragIter"""
     else:
       print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
-      return dragIter
+      return aeroCoeffs
 
 
   def RunPolar(self, alphaSeq):
diff --git a/dart/pyVII/vCoupler2.py b/vii/pyVII/vCoupler2.py
similarity index 100%
rename from dart/pyVII/vCoupler2.py
rename to vii/pyVII/vCoupler2.py
diff --git a/dart/pyVII/vInterpolator.py b/vii/pyVII/vInterpolator.py
similarity index 100%
rename from dart/pyVII/vInterpolator.py
rename to vii/pyVII/vInterpolator.py
diff --git a/dart/pyVII/vRbf.py b/vii/pyVII/vRbf.py
similarity index 100%
rename from dart/pyVII/vRbf.py
rename to vii/pyVII/vRbf.py
diff --git a/dart/pyVII/vUtils.py b/vii/pyVII/vUtils.py
similarity index 96%
rename from dart/pyVII/vUtils.py
rename to vii/pyVII/vUtils.py
index 7009f8d..16b8868 100644
--- a/dart/pyVII/vUtils.py
+++ b/vii/pyVII/vUtils.py
@@ -1,4 +1,5 @@
-import dart
+import vii
+from fwk.wutils import parseargs
 from fwk.coloring import ccolors
 
 def initBL(Re, Minf, CFL0, nSections, verbose=1):
@@ -16,8 +17,9 @@ def initBL(Re, Minf, CFL0, nSections, verbose=1):
     if not(0<=verbose<=3):
         print(ccolors.ANSI_RED, "dart::vUtils Fatal error : verbose not in valid range.", verbose, ccolors.ANSI_RESET)
         raise RuntimeError("Invalid parameter")
-    
-    solver = dart.ViscSolver(Re, Minf, CFL0, nSections, verbose)
+
+    args = parseargs()
+    solver = vii.ViscSolver(Re, Minf, CFL0, nSections, verbose=args.verb+1)
     return solver
 
 
diff --git a/vii/src/CMakeLists.txt b/vii/src/CMakeLists.txt
new file mode 100644
index 0000000..1e17fcd
--- /dev/null
+++ b/vii/src/CMakeLists.txt
@@ -0,0 +1,27 @@
+# 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.
+
+# CMake input file of dart.so
+
+FILE(GLOB SRCS *.h *.cpp *.inl *.hpp)
+
+ADD_LIBRARY(vii SHARED ${SRCS})
+MACRO_DebugPostfix(vii)
+TARGET_INCLUDE_DIRECTORIES(vii PUBLIC ${PROJECT_SOURCE_DIR}/vii/src)
+
+TARGET_LINK_LIBRARIES(vii tbox)
+
+INSTALL(TARGETS vii DESTINATION ${CMAKE_INSTALL_PREFIX})
+
+SOURCE_GROUP(base REGULAR_EXPRESSION ".*\\.(cpp|inl|hpp|h)")
diff --git a/vii/src/vii.h b/vii/src/vii.h
new file mode 100644
index 0000000..74bd31d
--- /dev/null
+++ b/vii/src/vii.h
@@ -0,0 +1,50 @@
+/*
+ * 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.
+ */
+
+// global header of the "vii" module
+
+#ifndef vii_H
+#define vii_H
+
+#if defined(WIN32)
+#ifdef vii_EXPORTS
+#define DART_API __declspec(dllexport)
+#else
+#define DART_API __declspec(dllimport)
+#endif
+#else
+#define DART_API
+#endif
+
+#include "tbox.h"
+
+/**
+ * @brief Namespace for dart module
+ */
+namespace vii
+{
+// Viscous
+class ViscSolver;
+class BLRegion;
+class TimeSolver;
+class ViscSolver;
+class Closures;
+class Interpolator;
+class ViscMesh;
+class InvBnd;
+} // namespace vii
+
+#endif //VII_H
diff --git a/dart/src/wBLRegion.cpp b/vii/src/wBLRegion.cpp
similarity index 99%
rename from dart/src/wBLRegion.cpp
rename to vii/src/wBLRegion.cpp
index f87aadb..0f83ecc 100644
--- a/dart/src/wBLRegion.cpp
+++ b/vii/src/wBLRegion.cpp
@@ -9,7 +9,7 @@
 #include <iomanip>
 #include <deque>
 
-using namespace dart;
+using namespace vii;
 
 BLRegion::BLRegion(double _xtrF)
 {
diff --git a/dart/src/wBLRegion.h b/vii/src/wBLRegion.h
similarity index 98%
rename from dart/src/wBLRegion.h
rename to vii/src/wBLRegion.h
index afd73d5..b109fdf 100644
--- a/dart/src/wBLRegion.h
+++ b/vii/src/wBLRegion.h
@@ -1,14 +1,14 @@
 #ifndef WBLREGION_H
 #define WBLREGION_H
 
-#include "dart.h"
+#include "vii.h"
 #include "wClosures.h"
 #include "wViscMesh.h"
 #include "wInvBnd.h"
 #include <vector>
 #include <string>
 
-namespace dart
+namespace vii
 {
 
 class DART_API BLRegion
diff --git a/dart/src/wClosures.cpp b/vii/src/wClosures.cpp
similarity index 99%
rename from dart/src/wClosures.cpp
rename to vii/src/wClosures.cpp
index 3e29976..10b2349 100644
--- a/dart/src/wClosures.cpp
+++ b/vii/src/wClosures.cpp
@@ -4,7 +4,7 @@
 #include <string>
 
 using namespace tbox;
-using namespace dart;
+using namespace vii;
 
 // This is an empty constructor
 Closures::Closures() {}
diff --git a/dart/src/wClosures.h b/vii/src/wClosures.h
similarity index 88%
rename from dart/src/wClosures.h
rename to vii/src/wClosures.h
index e86fc18..23f92d9 100644
--- a/dart/src/wClosures.h
+++ b/vii/src/wClosures.h
@@ -1,11 +1,11 @@
 #ifndef WCLOSURES_H
 #define WCLOSURES_H
 
-#include "dart.h"
+#include "vii.h"
 #include <vector>
 #include <string>
 
-namespace dart
+namespace vii
 {
 
 class DART_API Closures
@@ -16,5 +16,5 @@ public:
     std::vector<double> LaminarClosures(double theta, double H, double Ret, double Me, std::string name) const;
     std::vector<double> TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name) const;
 };
-} // namespace dart
+} // namespace vii
 #endif //WClosures_H
\ No newline at end of file
diff --git a/vii/src/wDartCoupler.cpp b/vii/src/wDartCoupler.cpp
new file mode 100644
index 0000000..cd40af4
--- /dev/null
+++ b/vii/src/wDartCoupler.cpp
@@ -0,0 +1,172 @@
+#include "wDartCoupler.h"
+#include "../../dart/src/wNewton.h"
+#include "../../dart/src/wBlowing.h"
+#include "wViscSolver.h"
+
+using namespace vii;
+
+DartCoupler::DartCoupler(const Newton *_invSolver, ViscSolver *_viscSolver, const std::vector<*Blowing> _blw, unsigned long _nCouplIter = 150, double _tol = 1e-4)
+{
+    invSolver = _invSolver;
+    viscSolver = _viscSolver;
+    blw = _blw;
+
+    maxCouplIter = _nCouplIter;
+    tol = _tol;
+
+    Initialize();
+}
+DartCoupler::~DartCoupler()
+
+void DartCoupler::Run(){
+    double cdPrev = 0.;
+    double error = 1.;
+
+    unsigned long nIters;
+
+    while (nIters < maxCouplIter){
+
+        invSolver.run();
+
+        ImposeInviscidBC();
+
+        viscSolver.Run();
+
+        error = (viscSolver.Cdt - cdPrev) / viscSolver.Cdt;
+
+        if (error < tol){
+            break;
+        }
+        else{
+            continue;
+        }
+
+        cdPrev = viscSolver.Cdt;
+
+        ImposeBlowing();
+
+        nIters++;
+
+    }
+
+    else{
+        std::cout << "Maximum number of " << maxCouplIter << " reached." << std::endl;
+    }
+}
+void DartCoupler::ImposeInviscidBC(){
+
+    for (int iReg=0; iReg < blw.size(); ++iReg){
+
+        for (auto iNode : blw[iReg]->nodes){
+            invSolver->U(iNode->row)
+        }
+    }
+}
+
+void DartCoupler::ImposeBlowing()
+
+void DartCoupler::Initialize(){
+
+    /* Create connectivity list */
+
+    std::vector<std::vector<unsigned int>,2> connectListNodes;
+    std::vector<std::vector<unsigned int>,2> connectListElems;
+
+    for (size_t iReg=0; iReg < blw.size(); ++iReg)
+    {
+        connectListNodes(iReg).resize(blw(iReg)->boundary->nodes.size(), 0);
+        connectListElems(iReg).resize(blw(iReg)->boudnary->tag->elems.size(), 0);
+    }
+
+
+    /* Find stagnation point */
+    int idxStag = 
+
+
+    N1 = np.zeros(self.nN, dtype=int) # Node number
+    connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
+    connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
+    data = np.zeros((self.boundary.nodes.size(), 10))
+    i = 0
+    for n in self.boundary.nodes:
+        data[i,0] = n.no
+        data[i,1] = n.pos[0]
+        data[i,2] = n.pos[1]
+        data[i,3] = n.pos[2]
+        data[i,4] = self.v[i,0]
+        data[i,5] = self.v[i,1]
+        data[i,6] = self.Me[i]
+        data[i,7] = self.rhoe[i]
+        data[i,8] = self.deltaStar[i]
+        data[i,9] = self.xx[i]
+        i += 1
+    # Table containing the element and its nodes
+    eData = np.zeros((self.nE,3), dtype=int)
+    for i in range(0, self.nE):
+        eData[i,0] = self.boundary.tag.elems[i].no
+        eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
+        eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
+    # Find the stagnation point
+    if couplIter == 0:
+        self.idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
+        idxStag = self.idxStag
+    else:
+        idxStag = self.idxStag
+    #idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
+    globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
+    # Find TE nodes (assuming suction side element is above pressure side element in the y direction)
+    idxTE = np.where(data[:,1] == np.max(data[:,1]))[0] # TE nodes
+    idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], eData[:,2] == data[idxTE[0],0]))[0] # TE element 1
+    idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], eData[:,2] == data[idxTE[1],0]))[0] # TE element 2
+    ye0 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe0[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe0[0],2])[0],2]) # y coordinates of TE element 1 CG
+    ye1 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe1[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe1[0],2])[0],2]) # y coordinates of TE element 2 CG
+    if ye0 - ye1 > 0: # element 1 containing TE node 1 is above element 2
+        upperGlobTE = int(data[idxTE[0],0])
+        lowerGlobTE = int(data[idxTE[1],0])
+    else:
+            upperGlobTE = int(data[idxTE[1],0])
+            lowerGlobTE = int(data[idxTE[0],0])
+    # Connectivity
+    connectListElems[0] = np.where(eData[:,1] == globStag)[0]
+    N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes -> globStag
+    connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
+    i = 1
+    upperTE = 0
+    lowerTE = 0
+    # Sort the suction part
+    while upperTE == 0:
+        N1[i] = eData[connectListElems[i-1],2] # Second node of the element
+        connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes
+        connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
+        if eData[connectListElems[i],2] == upperGlobTE:
+            upperTE = 1
+        i += 1
+    # Sort the pressure side
+    connectListElems[i] = np.where(eData[:,2] == globStag)[0]
+    connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
+    N1[i] = eData[connectListElems[i],2]
+    while lowerTE == 0:
+        N1[i+1]  = eData[connectListElems[i],1] # First node of the element
+        connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes
+        connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
+        if eData[connectListElems[i+1],1] == lowerGlobTE:
+            lowerTE = 1
+        i += 1
+    connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
+    data[:,0] = data[connectListNodes,0]
+    data[:,1] = data[connectListNodes,1]
+    data[:,2] = data[connectListNodes,2]
+    data[:,3] = data[connectListNodes,3]
+    data[:,4] = data[connectListNodes,4]
+    data[:,5] = data[connectListNodes,5]
+    data[:,6] = data[connectListNodes,6]
+    data[:,7] = data[connectListNodes,7]
+    data[:,8] = data[connectListNodes,8]
+    data[:,9] = data[connectListNodes,9]
+    # Separated upper and lower part
+    data = np.delete(data,0,1)
+    uData = data[0:np.argmax(data[:,0])+1]
+    lData = data[np.argmax(data[:,0])+1:None]
+    lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
+    return connectListNodes, connectListElems,[uData, lData]
+}
diff --git a/vii/src/wDartCoupler.h b/vii/src/wDartCoupler.h
new file mode 100644
index 0000000..8910ecd
--- /dev/null
+++ b/vii/src/wDartCoupler.h
@@ -0,0 +1,35 @@
+#ifndef WDARTCOUPLER_H
+#define WDARTCOUPLER_H
+
+#include "vii.h"
+#include "../../dart/src/wNewton.h"
+#include "../../dart/src/wBlowing.h"
+#include "wViscSolver.h"
+
+namespace vii
+{
+
+class DART_API DartCoupler
+{
+
+private:
+    unsigned long maxCouplIter;
+    double tol;
+
+
+public:
+
+    Newton *invSolver;
+    ViscSolver *viscSolver;
+    std::vector<Blowing> *blw;
+
+
+    DartCoupler();
+    ~DartCoupler();
+    Run();
+    ImposeBlowing();
+    ImposeInviscidBC();
+
+};
+} // namespace dart
+#endif //WDARTCOUPLER_H
diff --git a/dart/src/wInvBnd.cpp b/vii/src/wInvBnd.cpp
similarity index 93%
rename from dart/src/wInvBnd.cpp
rename to vii/src/wInvBnd.cpp
index d3ed14e..8d87227 100644
--- a/dart/src/wInvBnd.cpp
+++ b/vii/src/wInvBnd.cpp
@@ -21,7 +21,7 @@
 
 /* --- Namespaces -------------------------------------------------------------------- */
 
-using namespace dart;
+using namespace vii;
 
 /* ----------------------------------------------------------------------------------- */
 
@@ -33,7 +33,7 @@ double InvBnd::GetMaxMach() const
 {
     if (Me.size() <= 0)
     {
-        throw std::runtime_error("dart::InvBnd Invalid Mach number vector.");
+        throw std::runtime_error("vii::InvBnd Invalid Mach number vector.");
     }
     double Mmax = 0.0;
     for (size_t i=0; i<Me.size(); ++i)
diff --git a/dart/src/wInvBnd.h b/vii/src/wInvBnd.h
similarity index 95%
rename from dart/src/wInvBnd.h
rename to vii/src/wInvBnd.h
index d7f4299..b3f714c 100644
--- a/dart/src/wInvBnd.h
+++ b/vii/src/wInvBnd.h
@@ -1,12 +1,12 @@
 #ifndef WINVBND_H
 #define WINVBND_H
 
-#include "dart.h"
+#include "vii.h"
 #include "wBLRegion.h"
 #include <vector>
 #include <iostream>
 
-namespace dart
+namespace vii
 {
 
 class DART_API InvBnd
@@ -40,5 +40,5 @@ class DART_API InvBnd
     double GetMaxMach() const;
 
 };
-} // namespace dart
+} // namespace vii
 #endif //WINVBND_H
\ No newline at end of file
diff --git a/dart/src/wTimeSolver.cpp b/vii/src/wTimeSolver.cpp
similarity index 99%
rename from dart/src/wTimeSolver.cpp
rename to vii/src/wTimeSolver.cpp
index 3cae23f..dd42013 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/vii/src/wTimeSolver.cpp
@@ -9,7 +9,7 @@
 
 using namespace Eigen;
 using namespace tbox;
-using namespace dart;
+using namespace vii;
 
 TimeSolver::TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt, double _tol, unsigned int _itFrozenJac)
 {
diff --git a/dart/src/wTimeSolver.h b/vii/src/wTimeSolver.h
similarity index 97%
rename from dart/src/wTimeSolver.h
rename to vii/src/wTimeSolver.h
index 6ab52df..af485d4 100644
--- a/dart/src/wTimeSolver.h
+++ b/vii/src/wTimeSolver.h
@@ -1,12 +1,12 @@
 #ifndef WTIMESOLVER_H
 #define WTIMESOLVER_H
 
-#include "dart.h"
+#include "vii.h"
 #include "wViscFluxes.h"
 #include <vector>
 #include <string>
 
-namespace dart
+namespace vii
 {
 
 class DART_API TimeSolver
diff --git a/dart/src/wViscFluxes.cpp b/vii/src/wViscFluxes.cpp
similarity index 99%
rename from dart/src/wViscFluxes.cpp
rename to vii/src/wViscFluxes.cpp
index c940708..7fe7947 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/vii/src/wViscFluxes.cpp
@@ -9,7 +9,7 @@
 #include <algorithm>
 #include <string>
 
-using namespace dart;
+using namespace vii;
 using namespace Eigen;
 
 ViscFluxes::ViscFluxes(double _Re)
diff --git a/dart/src/wViscFluxes.h b/vii/src/wViscFluxes.h
similarity index 91%
rename from dart/src/wViscFluxes.h
rename to vii/src/wViscFluxes.h
index c295c98..a60bf35 100644
--- a/dart/src/wViscFluxes.h
+++ b/vii/src/wViscFluxes.h
@@ -1,13 +1,13 @@
 #ifndef WViscFluxes_H
 #define WViscFluxes_H
 
-#include "dart.h"
+#include "vii.h"
 #include "wBLRegion.h"
 #include <Eigen/Dense>
 #include <vector>
 #include <string>
 
-namespace dart
+namespace vii
 {
 
 class DART_API ViscFluxes
@@ -25,5 +25,5 @@ private:
     Eigen::VectorXd BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u);
     double AmplificationRoutine(double hk, double ret, double theta) const;
 };
-} // namespace dart
+} // namespace vii
 #endif //WViscFluxes_H
\ No newline at end of file
diff --git a/dart/src/wViscMesh.cpp b/vii/src/wViscMesh.cpp
similarity index 95%
rename from dart/src/wViscMesh.cpp
rename to vii/src/wViscMesh.cpp
index aafc154..9cb0544 100644
--- a/dart/src/wViscMesh.cpp
+++ b/vii/src/wViscMesh.cpp
@@ -17,7 +17,7 @@
 
 /* --- Namespaces -------------------------------------------------------------------- */
 
-using namespace dart;
+using namespace vii;
 
 /* ----------------------------------------------------------------------------------- */
 
@@ -29,7 +29,7 @@ void ViscMesh::SetGlob(std::vector<double> _x, std::vector<double> _y, std::vect
 {
   if (_x.size() != _y.size() || _y.size() != _z.size() || _x.size() != _z.size())
   {
-    throw std::runtime_error("dart::ViscMesh Wrong mesh sizes.");
+    throw std::runtime_error("vii::ViscMesh Wrong mesh sizes.");
   }
 
   nMarkers = _x.size();
diff --git a/dart/src/wViscMesh.h b/vii/src/wViscMesh.h
similarity index 97%
rename from dart/src/wViscMesh.h
rename to vii/src/wViscMesh.h
index ff1a53e..0531097 100644
--- a/dart/src/wViscMesh.h
+++ b/vii/src/wViscMesh.h
@@ -1,12 +1,12 @@
 #ifndef WVISCMESH_H
 #define WVISCMESH_H
 
-#include "dart.h"
+#include "vii.h"
 #include "wBLRegion.h"
 #include <vector>
 #include <iostream>
 
-namespace dart
+namespace vii
 {
 
 class DART_API ViscMesh
diff --git a/dart/src/wViscSolver.cpp b/vii/src/wViscSolver.cpp
similarity index 99%
rename from dart/src/wViscSolver.cpp
rename to vii/src/wViscSolver.cpp
index 1dc6e70..e4270fe 100644
--- a/dart/src/wViscSolver.cpp
+++ b/vii/src/wViscSolver.cpp
@@ -13,7 +13,6 @@
 #include "wViscSolver.h"
 #include "wBLRegion.h"
 #include "wTimeSolver.h"
-#include "wInterpolator.h"
 #include <iostream>
 #include <iomanip>
 #include <vector>
@@ -21,7 +20,7 @@
 
 /* --- Namespaces -------------------------------------------------------------------- */
 
-using namespace dart;
+using namespace vii;
 
 /* ----------------------------------------------------------------------------------- */
 
@@ -42,7 +41,6 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
 
     verbose = _verbose;
 
-
     /* Initialzing boundary layer */
 
     Sections.resize(nSections);
diff --git a/dart/src/wViscSolver.h b/vii/src/wViscSolver.h
similarity index 98%
rename from dart/src/wViscSolver.h
rename to vii/src/wViscSolver.h
index fdf8726..d500650 100644
--- a/dart/src/wViscSolver.h
+++ b/vii/src/wViscSolver.h
@@ -1,9 +1,9 @@
 #ifndef WVISCSOLVER_H
 #define WVISCSOLVER_H
 
-#include "dart.h"
+#include "vii.h"
 #include "wBLRegion.h"
-namespace dart
+namespace vii
 {
 
 class DART_API ViscSolver
diff --git a/dart/tests/bli_dualMesh.py b/vii/tests/bli2.py
similarity index 62%
rename from dart/tests/bli_dualMesh.py
rename to vii/tests/bli2.py
index ac96d71..c46a10a 100644
--- a/dart/tests/bli_dualMesh.py
+++ b/vii/tests/bli2.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python3 
 # -*- coding: utf-8 -*-
 
 # Copyright 2020 University of Liège
@@ -34,18 +34,18 @@
 # This test is provided to ensure that the solver works properly.
 # Mesh refinement may have to be performed to obtain physical results.
 
-import dart.utils as floU
-import dart.default as floD
-import dart.VII.viscUtils as viscU
+import dart.utils as invUtils
+import dart.default as invDefault
 
-import dart.VII.vCoupler as floVC
-import numpy as np
+import vii.pyVII.vCoupler as viiCoupler
+import vii.pyVII.vUtils as viscUtils
 
 import tbox
 import tbox.utils as tboxU
 import fwk
 from fwk.testing import *
 from fwk.coloring import ccolors
+from matplotlib import pyplot as plt
 
 import math
 
@@ -57,15 +57,11 @@ def main():
 
     # define flow variables
     Re = 1e7
-    alpha = 2*math.pi/180
-    #alphaSeq = np.arange(-5, 5.5, 0.5)
+    alpha = 5.*math.pi/180
     M_inf = 0.
 
-    meshFactor = 2
     CFL0 = 1
-
-    plotVar = 'H'
-
+    
     # ========================================================================================
 
     c_ref = 1
@@ -74,13 +70,14 @@ 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.001}
-    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+    msh, gmshWriter = invDefault.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
     # set the problem and add medium, initial/boundary conditions
     tms['pre'].start()
-    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+    pbl, _, _, bnd, blw = invDefault.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
     tms['pre'].stop()
 
     # solve viscous problem
@@ -88,53 +85,52 @@ def main():
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
     tms['solver'].start()
     
-    isolver = floD.newton(pbl)
-    vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
+    iSolver = invDefault.newton(pbl)
+    vSolver = viscUtils.initBL(Re, M_inf, CFL0, 1)
 
-    coupler = floVC.Coupler(isolver, vsolver)
-
-    coupler.CreateProblem(blw[0], blw[1])
-
-    #coupler.RunPolar(alphaSeq)
-    coupler.Run()
-    tms['solver'].stop()
+    Coupler = viiCoupler.Coupler(iSolver, vSolver)
+    Coupler.CreateProblem(blw[0], blw[1])
+    Coupler.Run()
 
     # extract Cp
-    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
-    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    Cp = invUtils.extract(bnd.groups[0].tag.elems, iSolver.Cp)
+    tboxU.write(Cp, 'Cp.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+
+    print('')
     # 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, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.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, iSolver.Cl, vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolver.Cm))
+
+    # Write results to file
+    vSolution = viscUtils.GetSolution(0, vSolver)
+    viscUtils.WriteFile(vSolution, vSolver.GetRe())
 
+    print('')
     # display timers
     tms['total'].stop()
     print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
     print('CPU statistics')
     print(tms)
-    
-    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)
-
-    blScalars, blVectors = viscU.ExtractVars(vsolver)
-    xtr=vsolver.Getxtr()
-    wData = viscU.Dictionary(blScalars, blVectors, xtr)
-    viscU.PlotVars(wData, plotVar)"""
+    print('SOLVERS statistics')
+    print(Coupler.tms)
 
     
     # check results
     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', isolver.Cl, 0.22113, 5e-3))
-        tests.add(CTest('Cd', vsolver.Cdt, 0.00534, 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))
+        tests.add(CTest('Cl', iSolver.Cl, 0.2208, 5e-3))
+        tests.add(CTest('Cd', vSolution['Cdt'], 0.00531, 1e-3, forceabs=True))
+        tests.add(CTest('Cdp', vSolution['Cdp'], 0.0009, 1e-3, forceabs=True))
+        tests.add(CTest('xtr_top', vSolution['xtrT'], 0.20, 0.03, forceabs=True))
+        tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.50, 0.03, forceabs=True))
+    if Re == 1e7 and M_inf == 0. and alpha == 5.*math.pi/180:
+        tests.add(CTest('Cl', iSolver.Cl, 0.554, 5e-3))
+        tests.add(CTest('Cd', vSolution['Cdt'], 0.007, 1e-3, forceabs=True))
+        tests.add(CTest('Cdp', vSolution['Cdp'],0.003, 1e-3, forceabs=True))
+        tests.add(CTest('xtr_top', vSolution['xtrT'], 0.06, 0.03, forceabs=True))
+        tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.74, 0.03, forceabs=True))
     else:
         raise Exception('Test not defined for this flow')
 
diff --git a/dart/tests/bli3.py b/vii/tests/bli3.py
similarity index 100%
rename from dart/tests/bli3.py
rename to vii/tests/bli3.py
-- 
GitLab