diff --git a/blast/interfaces/DDataStructure.py b/blast/interfaces/DDataStructure.py
index a5049f6ce9f6cef58c64e9b45f064ef73bc90264..93598c7c7008bc7cefa64894269f0e956e0c24fe 100644
--- a/blast/interfaces/DDataStructure.py
+++ b/blast/interfaces/DDataStructure.py
@@ -14,7 +14,7 @@ class Group:
         self.nDim = _nDim
         self.stagPoint = None
 
-    def initStructures(self, nPoints):
+    def initStructures(self, nPoints, nElems):
         """Initialize the data structures.
         
         args:
@@ -23,7 +23,7 @@ class Group:
                 Number of points in the region
         """
         self.nPoints = nPoints
-        self.nElems = nPoints - 1
+        self.nElems = nElems
         self.nodesCoord = np.zeros((self.nPoints,4))
         
         self.V = np.zeros((self.nPoints,3))
diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index 29e4081378ddbe88b673f6f82a8d1e312da5df70..a5aa22b140d9ebe950ca197c2c3dfad04142191d 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/blast/interfaces/DSolversInterface.py
@@ -22,7 +22,7 @@ class SolversInterface:
             # Initialize viscous structures for matching computations
             for iSec in range(len(self.vBnd)):
                 for iReg, reg in enumerate(self.vBnd[iSec]):
-                    reg.initStructures(self.iBnd[iReg].nPoints)
+                    reg.initStructures(self.iBnd[iReg].nPoints, self.iBnd[iReg].nElems)
                     self.vBnd[iSec][iReg].nodesCoord = self.iBnd[iReg].nodesCoord.copy()
                     self.vBnd[iSec][iReg].elemsCoord = self.iBnd[iReg].elemsCoord.copy()
                     self.vBnd[iSec][iReg].connectListNodes = self.iBnd[iReg].connectListNodes.copy()
@@ -306,7 +306,8 @@ class SolversInterface:
     def getVWing2(self):
         from scipy.interpolate import griddata
         from scipy.interpolate import interp1d
-        import matplotlib.pyplot as plt
+
+        ndim = self.getnDim()
 
         upper_elems = []
         lower_elems = []
@@ -320,20 +321,23 @@ class SolversInterface:
         for e in self.blw[0].tag.elems:
             if e.no in lower_elems:
                 for n in e.nodes:
-                    if n.row not in lower_nodes[:,3]:
+                    if n.row not in lower_nodes[:,-1]:
                         lower_nodes = np.row_stack((lower_nodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
             else:
                 upper_elems.append(e.no)
                 for n in e.nodes:
-                    if n.row not in upper_nodes[:,3]:
+                    if n.row not in upper_nodes[:,-1]:
                         upper_nodes = np.row_stack((upper_nodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
 
+        print('Found upper side with', len(upper_elems), 'elements')
+        print('Found lower side with', len(lower_elems), 'elements')
         print('Found upper side with', len(upper_nodes), 'points')
         print('Found lower side with', len(lower_nodes), 'points')
         
         # Remove points higher than span
-        upper_nodes = upper_nodes[upper_nodes[:,1] < self.cfg['span']*0.99, :]
-        lower_nodes = lower_nodes[lower_nodes[:,1] < self.cfg['span']*0.99, :]
+        if ndim == 3:
+            upper_nodes = upper_nodes[upper_nodes[:,1] < self.cfg['span']*0.99, :]
+            lower_nodes = lower_nodes[lower_nodes[:,1] < self.cfg['span']*0.99, :]
 
         if upper_nodes.shape[0] == 0 or lower_nodes.shape[0] == 0:
             raise RuntimeError('Could not identify upper or lower side.')
@@ -343,46 +347,41 @@ class SolversInterface:
         interpolated_lower = []
 
         for y in sections:
-            x = np.linspace(np.min(upper_nodes[:, 0]), np.max(upper_nodes[:, 0]), 1000)
-            z_upper = griddata(upper_nodes[:, :2], upper_nodes[:, 2], (x, np.full_like(x, y)), method='cubic')
-            z_lower = griddata(lower_nodes[:, :2], lower_nodes[:, 2], (x, np.full_like(x, y)), method='cubic')
-            interpolated_upper.append(np.column_stack((x, np.full_like(x, y), z_upper)))
-            interpolated_lower.append(np.column_stack((x, np.full_like(x, y), z_lower)))
-            interpolated_upper[-1] = np.delete(interpolated_upper[-1], np.where(np.isnan(interpolated_upper[-1][:, 2])), axis=0)
-            interpolated_lower[-1] = np.delete(interpolated_lower[-1], np.where(np.isnan(interpolated_lower[-1][:, 2])), axis=0)
-
-        """# Ensure that we have one LE point after interpolation
-        for i in range(len(interpolated_upper)):
-            avrg_le = (interpolated_upper[i][0,2] + interpolated_lower[i][0,2]) / 2
-            interpolated_upper[i][0,2] = avrg_le
-            interpolated_lower[i][0,2] = avrg_le
-            avrg_te = (interpolated_upper[i][-1,2] + interpolated_lower[i][-1,2]) / 2
-            interpolated_upper[i][-1,2] = avrg_te
-            interpolated_lower[i][-1,2] = avrg_te"""
-        
+            data_up = np.zeros((self.cfg['nPoints'], 3))
+            data_lw = np.zeros((self.cfg['nPoints'], 3))
+            data_up[:,0] = np.linspace(np.min(upper_nodes[:, 0]), np.max(upper_nodes[:, 0]), self.cfg['nPoints'])
+            data_lw[:,0] = np.linspace(np.min(upper_nodes[:, 0]), np.max(upper_nodes[:, 0]), self.cfg['nPoints'])
+            data_up[:,1] = y * np.ones(data_up.shape[0])
+            data_lw[:,1] = y * np.ones(data_lw.shape[0])
+            if ndim == 2:
+                data_up[:,ndim-1] = interp1d(upper_nodes[:, 0], upper_nodes[:, ndim-1], kind='cubic')(data_up[:,0])
+                data_lw[:,ndim-1] = interp1d(lower_nodes[:, 0], lower_nodes[:, ndim-1], kind='cubic')(data_lw[:,0])
+            else:
+                data_up[:,ndim-1] = griddata(upper_nodes[:, :ndim-1], upper_nodes[:, ndim-1], (data_up[:,0], np.full_like(data_up[:,0], y)), method='cubic')
+                data_lw[:,ndim-1] = griddata(lower_nodes[:, :ndim-1], lower_nodes[:, ndim-1], (data_lw[:,0], np.full_like(data_lw[:,0], y)), method='cubic')
+            interpolated_upper.append(data_up)
+            interpolated_lower.append(data_lw)
+            interpolated_upper[-1] = np.delete(interpolated_upper[-1], np.where(np.isnan(interpolated_upper[-1][:, ndim-1])), axis=0)
+            interpolated_lower[-1] = np.delete(interpolated_lower[-1], np.where(np.isnan(interpolated_lower[-1][:, ndim-1])), axis=0)
+
         # Create sections
         sections_final = []
         zeta = np.linspace(0., np.pi, self.cfg['nPoints'])
         for i in range(len(interpolated_upper)):
             chord_upper = np.max(interpolated_upper[i][:,0]) - np.min(interpolated_upper[i][:,0])
-            xu = ( chord_upper / 2) * (np.cos(zeta) + 1) + np.min(interpolated_upper[i][:,0])
-            zu = interp1d(interpolated_upper[i][:,0], interpolated_upper[i][:,2], kind='cubic', fill_value='extrapolate')(xu)
+            data_up = np.zeros((self.cfg['nPoints'], 3))
+            data_up[:,0] = ( chord_upper / 2) * (np.cos(zeta) + 1) + np.min(interpolated_upper[i][:,0])
+            data_up[:,1] = sections[i] * np.ones(data_up.shape[0])
+            data_up[:,ndim-1] = interp1d(interpolated_upper[i][:,0], interpolated_upper[i][:,ndim-1], kind='cubic', fill_value='extrapolate')(data_up[:,0])
 
             chord_lower = np.max(interpolated_lower[i][:,0]) - np.min(interpolated_lower[i][:,0])
-            xl = (chord_lower / 2) * (np.cos(zeta) + 1) + np.min(interpolated_lower[i][:,0])
-            zl = interp1d(interpolated_lower[i][:,0], interpolated_lower[i][:,2], kind='cubic', fill_value='extrapolate')(xl)
 
-            uf = np.column_stack((xu, np.full_like(xu, sections[i]), zu))
-            lf = np.column_stack((xl, np.full_like(xl, sections[i]), zl))
+            data_lw = np.zeros((self.cfg['nPoints'], 3))
+            data_lw[:,0] = (chord_lower / 2) * (np.cos(zeta) + 1) + np.min(interpolated_lower[i][:,0])
+            data_lw[:,1] = sections[i] * np.ones(data_lw.shape[0])
+            data_lw[:,ndim-1] = interp1d(interpolated_lower[i][:,0], interpolated_lower[i][:,ndim-1], kind='cubic', fill_value='extrapolate')(data_lw[:,0])
 
-            sections_final.append(np.row_stack((uf, np.flip(lf, axis=0)[1:])))
-
-        # Compute elements positions
-        elems_coords_wing = []
-        for isec, sec in enumerate(sections_final):
-            elems_coords_wing.append(np.zeros((sec.shape[0]-1, 3)))
-            for i in range(sec.shape[0]-1):
-                elems_coords_wing[isec][i, :] = (sec[i, :] + sec[i+1, :]) / 2
+            sections_final.append(np.row_stack((data_up, np.flip(data_lw, axis=0)[1:])))
         
         # Wake
         wake = np.zeros((0,3))
@@ -393,29 +392,38 @@ class SolversInterface:
         interpolated_wake = []
 
         for y in sections:
-            x = np.linspace(np.min(wake[:, 0]), np.max(wake[:, 0]), 1000)
-            z = griddata(wake[:, :2], wake[:, 2], (x, np.full_like(x, y)), method='cubic')
-            interpolated_wake.append(np.column_stack((x, np.full_like(x, y), z)))
-            interpolated_wake[-1] = np.delete(interpolated_wake[-1], np.where(np.isnan(interpolated_wake[-1][:, 2])), axis=0)
-
-        # Ensure that we have one LE point after interpolation
-        for i in range(len(interpolated_wake)):
-            avrg_le = interpolated_wake[i][0,2]
-            interpolated_wake[i][0,2] = avrg_le
-            avrg_te = interpolated_wake[i][-1,2]
-            interpolated_wake[i][-1,2] = avrg_te
+            data = np.zeros((self.cfg['nPoints'], 3))
+            data[:,0] = np.linspace(np.min(wake[:, 0]), np.max(wake[:, 0]), self.cfg['nPoints'])
+            # set all values in data[:,1] to y
+            data[:,1] = y * np.ones(data.shape[0])
+            if ndim == 2:
+                data[:, ndim-1] = interp1d(wake[:, 0], wake[:, ndim-1], kind='cubic')(data[:,0])
+            else:
+                data[:, ndim-1] = griddata(wake[:, :ndim-1], wake[:, ndim-1], (data[:,0], np.full_like(data[:,0], data[:,1])), method='cubic')
+            interpolated_wake.append(data)
+            interpolated_wake[-1] = np.delete(interpolated_wake[-1], np.where(np.isnan(interpolated_wake[-1][:, ndim-1])), axis=0)
 
         # Create sections
         wake_final = []
-
         for i in range(len(interpolated_wake)):
             chord = np.max(interpolated_wake[i][:,0]) - np.min(interpolated_wake[i][:,0])
             x = (chord / 2) * (np.cos(zeta) + 1) + np.min(interpolated_wake[i][:,0])
-            z = interp1d(interpolated_wake[i][:,0], interpolated_wake[i][:,2], kind='cubic')(x)
-
+            z = interp1d(interpolated_wake[i][:,0], interpolated_wake[i][:,ndim-1], kind='cubic')(x)
             wf = np.column_stack((x, np.full_like(x, sections[i]), z))
             wake_final.append(np.flip(wf, axis=0))
         
+        # Ensure that there is always only one TE/first wake point.
+        for i in range(len(sections_final)):
+            sections_final[i][0, :] = wake_final[i][0, :]
+            sections_final[i][-1, :] = wake_final[i][0, :]
+
+        # Compute elements positions
+        elems_coords_wing = []
+        for isec, sec in enumerate(sections_final):
+            elems_coords_wing.append(np.zeros((sec.shape[0]-1, sec.shape[1])))
+            for i in range(sec.shape[0]-1):
+                elems_coords_wing[isec][i, :] = (sec[i, :] + sec[i+1, :]) / 2
+        
         elems_coords_wake = []
         for isec, sec in enumerate(wake_final):
             elems_coords_wake.append(np.zeros((sec.shape[0]-1, 3)))
@@ -423,134 +431,16 @@ class SolversInterface:
                 elems_coords_wake[isec][i, :] = (sec[i, :] + sec[i+1, :]) / 2
 
         for iy in range(len(sections_final)):
-            self.vBnd[iy][0].initStructures(sections_final[iy].shape[0])
+            self.vBnd[iy][0].initStructures(sections_final[iy].shape[0], elems_coords_wing[iy].shape[0])
             self.vBnd[iy][0].nodesCoord = sections_final[iy]
             self.vBnd[iy][0].elemsCoord = elems_coords_wing[iy]
-
         for iy in range(len(wake_final)):
-            self.vBnd[iy][1].initStructures(wake_final[iy].shape[0])
+            self.vBnd[iy][1].initStructures(wake_final[iy].shape[0], elems_coords_wake[iy].shape[0])
             self.vBnd[iy][1].nodesCoord = wake_final[iy]
             self.vBnd[iy][1].elemsCoord = elems_coords_wake[iy]
-        
-        print('Created', len(sections_final), 'sections and', len(wake_final), 'wake sections')
-        # plot
-        fig, ax = plt.subplots()
-        ax = fig.add_subplot(111, projection='3d')
-        ax.plot_trisurf(upper_nodes[:,0], upper_nodes[:,1], upper_nodes[:,2], color='blue', alpha=0.5)
-        ax.plot_trisurf(lower_nodes[:,0], lower_nodes[:,1], lower_nodes[:,2], color='blue', alpha=0.5)
-        #ax.plot_trisurf(wake[:,0], wake[:,1], wake[:,2], color='green', alpha=0.5)
-        for iy in range(len(sections_final)):
-            ax.plot(sections_final[iy][:,0], sections_final[iy][:,1], sections_final[iy][:,2], '-')
-        for iy in range(len(wake_final)):
-            ax.plot(wake_final[iy][:,0], wake_final[iy][:,1], wake_final[iy][:,2], '-', color='red')
-            #ax.plot_trisurf(ib.nodesCoord[ib.nodesCoord[:,1]<0.99*self.cfg['span'],0], ib.nodesCoord[ib.nodesCoord[:,1]<0.99*self.cfg['span'],1], ib.nodesCoord[ib.nodesCoord[:,1]<0.99*self.cfg['span'],2], color='red', alpha=0.5)
-        # for sec in sections_final:
-        #     ax.plot(sec[:,0], sec[:,2])
-        # for sec in wake_final:
-        #     ax.plot(sec[:,0], sec[:,2])
-        plt.axis('equal')
-        plt.savefig("ax.png")
-        # plt.draw()
-        # plt.show()
-        # quit()
-
-        for i, sec in enumerate(self.vBnd):
-            plt.figure()
-            for reg in sec:
-                if reg.name == 'wake':
-                    plt.plot(reg.nodesCoord[:,0], reg.nodesCoord[:,2], '-')
-                else:
-                    plt.plot(reg.nodesCoord[:,0], reg.nodesCoord[:,2], '-')
-            plt.axis('equal')
-            plt.xlim([0, 1.2])
-            #plt.draw()
-            plt.savefig(f'Sec{i}.png')
-
-        return np.array(upper_nodes[:,3], dtype=int), np.array(lower_nodes[:,3], dtype=int), upper_elems, lower_elems
-
-    def getVWing(self):
-        """Obtain the nodes' location and row and cg of all elements.
-        If the mesh is the viscous one, TE nodes are duplicated (by creating a new row).
-        They exist for the interpolation and for the solution computation.
-        If not, the nodes are pseudo-duplicated (without changing the row). They therefore
-        exist only for the interpolation but not for the solution.
-        """
-        data = [np.zeros((0,4)) for _ in range(2)]
-        lwRows = []
-        for e in self.cfg['vMsh'].elems:
-            if e.ptag.name == 'wing' or e.ptag.name == 'wing_'\
-                or e.ptag.name == 'leadingEdge':
-                iRegion = 0
-            elif e.ptag.name == 'wake' or e.ptag.name == 'wake_':
-                iRegion = 1
-            else:
-                continue
-            
-            # Get nodes
-            for nw in e.nodes:
-                if nw.row not in data[iRegion][:,3]:
-                    data[iRegion] = np.vstack((data[iRegion], [nw.pos[0],\
-                                                                nw.pos[1],\
-                                                                nw.pos[2],\
-                                                                nw.row]))
-                    # If wing_ exists, the lower side is identifiable.
-                    if e.ptag.name == 'wing_':
-                        lwRows.append(nw.row)
-        
-        if len(lwRows) == 0:
-            raise RuntimeError('Could not identify lower side.')
-
-        self.vlwIdx = []
-        for idx in range(len(data[0][:,3])):
-            if data[0][idx, 3] in lwRows:
-                self.vlwIdx.append(idx)
-
-        maxRowNb = max(data[0][:,3])
-        # Initialize viscous strutures
-        self.cfg['EffSections'] = np.empty(0)
-
-        for iReg, reg in enumerate(data):
-            for iy, y in enumerate(self.cfg['sections']):
-                # Find nearest value in the mesh
-                try:
-                    idx = (np.abs(reg[:,1]-y).argmin())
-                except:
-                    print('Region', reg[:,1], 'cannot be found through\
-                          config input', y)
-                    raise RuntimeError("Could not identify section.\
-                                       Possibly a missing tag")
-                section = reg[idx,1]
-                if iReg == 0:
-                    self.cfg['EffSections'] = np.append(self.cfg['EffSections'], section)
-                nodesSec = reg[abs(reg[:,1]-section)<1e-3, :]
-
-                # Separate points on upper and lower side
-                if iReg == 0:
-                    lower_points = []
-                    upper_points = []
-                    for point in nodesSec:
-                        if point[3] in lwRows:
-                            lower_points.append(point)
-                        else:
-                            upper_points.append(point)
-
-                    # Sort points in selig format
-                    lower_points = sorted(lower_points, key=lambda p: p[0])
-                    upper_points = sorted(upper_points, key=lambda p: p[0], reverse=True)
-                    nodesSec = np.row_stack((upper_points, lower_points))
-                    nodesSec = np.row_stack((nodesSec, [nodesSec[0,0], nodesSec[0,1], nodesSec[0,2], nodesSec[0,3]+maxRowNb]))
-                else:
-                    nodesSec = nodesSec[nodesSec[:,0].argsort()]
 
-                # Compute cg pos on the section (!= cg of 2D elems on the surface)
-                cgSec = np.zeros((len(nodesSec[:,0])-1, 4))
-                for i in range(len(cgSec)):
-                    cgSec[i,:3] = (nodesSec[i,:3] + nodesSec[i+1,:3])/2
-                    cgSec[i, 3] = i
-
-                self.vBnd[iy][iReg].initStructures(nodesSec.shape[0])
-                self.vBnd[iy][iReg].nodesCoord = nodesSec
-                self.vBnd[iy][iReg].elemsCoord = cgSec
+        print('Created', len(sections_final), 'sections and', len(wake_final), 'wake sections')
+        return np.array(upper_nodes[:,3], dtype=int), np.array(lower_nodes[:,3], dtype=int), np.array(upper_elems, dtype=int), np.array(lower_elems, dtype=int)
 
     # Inviscid solver methods
     def run()->int:
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index ffa71953eba16204c3d6d41525f4603c1c7dd4b9..0fccfc1fce37b4f0269c2bde3f10ba7be26914e1 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/blast/interfaces/dart/DartInterface.py
@@ -179,7 +179,7 @@ class DartInterface(SolversInterface):
     def __getWing2D(self):
         """Create data structure for information transfer.
         """
-        self.iBnd[0].initStructures(self.blw[0].nodes.size())
+        self.iBnd[0].initStructures(self.blw[0].nodes.size(), self.blw[0].tag.elems.size())
         # Node number
         N1 = np.zeros(self.blw[0].nodes.size(), dtype=int)
         # Index in boundary.nodes
@@ -291,11 +291,11 @@ class DartInterface(SolversInterface):
             elemCoord[i,0] = 0.5 * (self.iBnd[0].nodesCoord[i,0] + self.iBnd[0].nodesCoord[i+1,0])
             elemCoord[i,1] = 0.5 * (self.iBnd[0].nodesCoord[i,1] + self.iBnd[0].nodesCoord[i+1,1])
             elemCoord[i,2] = 0.5 * (self.iBnd[0].nodesCoord[i,2] + self.iBnd[0].nodesCoord[i+1,2])
-            elemCoord[i,3] = i
+            elemCoord[i,3] = self.blw[0].tag.elems[i].no
         self.iBnd[0].elemsCoord = elemCoord
 
         # Wake
-        self.iBnd[1].initStructures(self.blw[1].nodes.size())
+        self.iBnd[1].initStructures(self.blw[1].nodes.size(), self.blw[1].tag.elems.size())
         # Node number
         N1 = np.zeros(self.blw[1].nodes.size(), dtype=int)
         # Index in boundary.nodes
@@ -376,18 +376,8 @@ class DartInterface(SolversInterface):
                 cg_pt/=e.nodes.size()
                 cg_pt = np.hstack((cg_pt, e.no))
                 cg[iRegion] = np.vstack((cg[iRegion], cg_pt))
-        # self.ilwIdx = []
-        # if 'Sym' in self.cfg:
-        #     for s in self.cfg['Sym']:
-        #         for iReg in range(len(data)):
-        #             dummy = data[iReg].copy()
-        #             dummy[:,1] -= s
-        #             dummy[:,1] *= -1
-        #             dummy[:,1] += s
-        #             data[iReg] = np.row_stack((data[iReg], dummy))
-
         for n in range(len(data)):
-            self.iBnd[n].initStructures(data[n].shape[0])
+            self.iBnd[n].initStructures(data[n].shape[0], cg[n].shape[0])
             self.iBnd[n].nodesCoord = data[n]
             self.iBnd[n].elemsCoord = cg[n]
             self.iBnd[n].setConnectList(data[n][:,3], cg[n][:,3])
diff --git a/blast/interfaces/interpolators/DRbfInterpolator.py b/blast/interfaces/interpolators/DRbfInterpolator.py
index 4f5ae028116fb80d3fb6f0591660de1c0fe20816..e27ababc1ae40fd46962194d692f31a53cd8e7af 100644
--- a/blast/interfaces/interpolators/DRbfInterpolator.py
+++ b/blast/interfaces/interpolators/DRbfInterpolator.py
@@ -15,10 +15,8 @@ class RbfInterpolator:
 
 
     def inviscidToViscous(self, iDict, vDict):
-        
-        print('Inviscid to viscous')
         # upper side
-        x0_up = np.zeros((0, self.nDim+1))
+        x0_up = np.zeros((0, 4))
         for row in self.upNodes:
             x0_up = np.row_stack((x0_up, np.unique(iDict[0].nodesCoord[iDict[0].nodesCoord[:,3]==row,:], axis=0)))
         M0_up = np.zeros(x0_up.shape[0])
@@ -32,7 +30,7 @@ class RbfInterpolator:
             for idim in range(self.nDim):
                 v0_up[i, idim] = iDict[0].V[iDict[0].nodesCoord[:,3] == row, idim][0]
         # lower side
-        x0_lw = np.zeros((0, self.nDim+1))
+        x0_lw = np.zeros((0, 4))
         for row in self.lwNodes:
             x0_lw = np.row_stack((x0_lw, np.unique(iDict[0].nodesCoord[iDict[0].nodesCoord[:,3]==row,:], axis=0)))
         M0_lw = np.zeros(x0_lw.shape[0])
@@ -47,7 +45,6 @@ class RbfInterpolator:
                 v0_lw[i, idim] = iDict[0].V[iDict[0].nodesCoord[:,3] == row, idim][0]
 
         for iSec in range(len(vDict)):
-
             xInterp_up = vDict[iSec][0].nodesCoord[:np.argmin(vDict[iSec][0].nodesCoord[:,0]), :self.nDim]
             M_up = self.__rbfToSection(x0_up[:,:self.nDim], M0_up, xInterp_up)
             rho_up = self.__rbfToSection(x0_up[:,:self.nDim], rho0_up, xInterp_up)
@@ -76,7 +73,35 @@ class RbfInterpolator:
             vDict[iSec][0].updateVars(v_int, M_int, rho_int)
 
             vDict[iSec][1].updateVars(v_wk, M_wk, rho_wk)
-            
+
+        # from matplotlib import pyplot as plt
+        # for iSec in range(len(vDict)):
+        #     plt.figure()
+        #     plt.plot(iDict[0].nodesCoord[:,0], iDict[0].M, '-', label='inv')
+        #     plt.plot(iDict[1].nodesCoord[:,0], iDict[1].M, '-', label='inv wk')
+        #     plt.plot(vDict[iSec][0].nodesCoord[:,0], vDict[iSec][0].M, '--', label='visc')
+        #     plt.plot(vDict[iSec][1].nodesCoord[:,0], vDict[iSec][1].M, '--', label='visc wk')
+        #     plt.legend(frameon=False)
+        #     plt.title('M')
+        #     plt.draw()
+        #     plt.figure()
+        #     plt.plot(iDict[0].nodesCoord[:,0], iDict[0].Rho, '-', label='inv')
+        #     plt.plot(iDict[1].nodesCoord[:,0], iDict[1].Rho, '-', label='inv wk')
+        #     plt.plot(vDict[iSec][0].nodesCoord[:,0], vDict[iSec][0].Rho, '--', label='visc') 
+        #     plt.plot(vDict[iSec][1].nodesCoord[:,0], vDict[iSec][1].Rho, '--', label='visc wk')
+        #     plt.legend(frameon=False)
+        #     plt.title('Rho')
+        #     plt.draw()
+        #     for idim in range(self.nDim):
+        #         plt.figure()
+        #         plt.plot(iDict[0].nodesCoord[:,0], iDict[0].V[:,idim], '-', label='inv')
+        #         plt.plot(iDict[1].nodesCoord[:,0], iDict[1].V[:,idim], '-', label='inv wk')
+        #         plt.plot(vDict[iSec][0].nodesCoord[:,0], vDict[iSec][0].V[:,idim], '--', label='visc')
+        #         plt.plot(vDict[iSec][1].nodesCoord[:,0], vDict[iSec][1].V[:,idim], '--', label='visc wk')
+        #         plt.legend(frameon=False)
+        #         plt.title('V' + str(idim))
+        #         plt.draw()
+        # plt.show()
 
             # M_tot = self.__rbfToSection(iDict[0].nodesCoord[:,:self.nDim], iDict[0].M, vDict[iSec][0].nodesCoord[:,:self.nDim])
             # V_tot = np.zeros((vDict[iSec][0].nodesCoord.shape[0], self.nDim))
@@ -122,37 +147,8 @@ class RbfInterpolator:
             #     plt.title('Velocity'+str(idim))
             #     plt.draw()
             # plt.show()
-            
-
-
-
-            # # Plot 3D x0_lw
-            # fig = plt.figure()
-            # ax = fig.add_subplot(111, projection='3d')
-            # ax.scatter(x0_lw[:,0], x0_lw[:,1], x0_lw[:,2])
-            # plt.show()
-
-
-
-            
-
-
-    def inviscidToViscous2(self, iDict, vDict):
-        ## Airfoil
-        # Find stagnation point
-        for iSec in range(len(vDict)):
-            for iReg, reg in enumerate(vDict[iSec]):
-                v = np.zeros((reg.nodesCoord.shape[0], 3))
-                M = np.zeros(reg.nodesCoord.shape[0])
-                rho = np.zeros(reg.nodesCoord.shape[0])
-                for iDim in range(3):
-                    v[:,iDim] = self.__rbfToSection(iDict[iReg].nodesCoord[:,:(self.cfg['nDim'] if iDict[iReg].name == 'iWing' else 1)], iDict[iReg].V[:,iDim], reg.nodesCoord[:,:(self.cfg['nDim'] if 'vAirfoil' in reg.name else 1)])
-                M = self.__rbfToSection(iDict[iReg].nodesCoord[:,:(self.cfg['nDim'] if iDict[iReg].name == 'iWing' else 1)], iDict[iReg].M, reg.nodesCoord[:,:(self.cfg['nDim'] if 'vAirfoil' in reg.name else 1)])
-                rho = self.__rbfToSection(iDict[iReg].nodesCoord[:,:(self.cfg['nDim'] if iDict[iReg].name == 'iWing' else 1)], iDict[iReg].Rho, reg.nodesCoord[:,:(self.cfg['nDim'] if 'vAirfoil' in reg.name else 1)])
-                vDict[iSec][iReg].updateVars(v, M, rho)
 
     def viscousToInviscid(self, iDict, vDict):
-
         # Viscous side
         xv_up = np.zeros((0, self.nDim))
         xv_lw = np.zeros((0, self.nDim))
@@ -169,60 +165,44 @@ class RbfInterpolator:
             blwv_wk = np.concatenate((blwv_wk, sec[1].blowingVel))
         
         # Inviscid side
-        xi_up = np.zeros((0, self.nDim+1))
-        xi_lw = np.zeros((0, self.nDim+1))
+        xi_up = np.zeros((0, 4))
+        xi_lw = np.zeros((0, 4))
         for no in self.upElems:
-            xi_up = np.row_stack((xi_up, np.unique(iDict[0].elemsCoord[iDict[0].elemsCoord[:,3]==no,:], axis=0)))
+            xi_up = np.row_stack((xi_up, np.unique(iDict[0].elemsCoord[iDict[0].elemsCoord[:,-1]==no,:], axis=0)))
         for no in self.lwElems:
-            xi_lw = np.row_stack((xi_lw, np.unique(iDict[0].elemsCoord[iDict[0].elemsCoord[:,3]==no,:], axis=0)))
+            xi_lw = np.row_stack((xi_lw, np.unique(iDict[0].elemsCoord[iDict[0].elemsCoord[:,-1]==no,:], axis=0)))
         xi_wk = iDict[1].elemsCoord[:,:self.nDim]
 
         blwi_up = self.__rbfToSection(xv_up, blwv_up, xi_up[:,:self.nDim])
         blwi_lw = self.__rbfToSection(xv_lw, blwv_lw, xi_lw[:,:self.nDim])
         blwi_wk = self.__rbfToSection(xv_wk, blwv_wk, xi_wk[:,:self.nDim])
 
-        iDict[0].blowingVel = np.concatenate((blwi_up, blwi_lw))
+        for no in xi_up[:,-1]:
+            idx = np.where(iDict[0].elemsCoord[:,3]==no)[0]
+            iDict[0].blowingVel[idx] = blwi_up[xi_up[:,-1]==no]
+        for no in xi_lw[:,-1]:
+            idx = np.where(iDict[0].elemsCoord[:,3]==no)[0]
+            iDict[0].blowingVel[idx] = blwi_lw[xi_lw[:,-1]==no]
         iDict[1].blowingVel = blwi_wk
 
-        """from matplotlib import pyplot as plt
-        from mpl_toolkits.mplot3d import Axes3D
-        # plot 3D x0_up
-        fig = plt.figure()
-        ax = fig.add_subplot(111, projection='3d')
-
-        ax.plot_trisurf(xi_up[:,0], xi_up[:,1], blwi_up, color='blue')
-        ax.plot(xv_up[:,0], xv_up[:,1], blwv_up, color='red')
-        plt.draw()
-
-        fig = plt.figure()
-        ax = fig.add_subplot(111, projection='3d')
-        ax.plot_trisurf(xi_lw[:,0], xi_lw[:,1], blwi_lw, color='blue')
-        ax.plot(xv_lw[:,0], xv_lw[:,1], blwv_lw, color='red')
-        plt.draw()
-
-        plt.show()"""
-
-    def viscousToInviscid2(self, iDict, vDict):
-        print('Viscous to inviscid')
-        if self.nDim == 2:
-            for iReg, reg in enumerate(iDict):
-                iDict[iReg].blowingVel = self.__rbfToSection(vDict[0][iReg].elemsCoordTr[:,:(self.cfg['nDim']-1 if 'vAirfoil' in reg.name else 1)], vDict[0][iReg].blowingVel, reg.elemsCoordTr[:,:(self.cfg['nDim']-1 if reg.name == 'iWing' else 1)])
-        elif self.nDim == 3:
-            for iReg, reg in enumerate(iDict):
-                viscElemsCoord = np.zeros((0,3))
-                viscBlowing = np.zeros(0)
-                for iSec, sec in enumerate(vDict):
-                    viscElemsCoord = np.row_stack((viscElemsCoord, sec[iReg].elemsCoord[:,:3]))
-                    viscBlowing = np.concatenate((viscBlowing, sec[iReg].blowingVel))
-                if 'Sym' in self.cfg:
-                    for s in self.cfg['Sym']:
-                        dummy = viscElemsCoord.copy()
-                        dummy[:,1] = (dummy[:,1] - s)*(-1) + s
-                        viscElemsCoord = np.row_stack((viscElemsCoord, dummy))
-                        viscBlowing = np.concatenate((viscBlowing, viscBlowing))
-                reg.blowingVel = self.__rbfToSection(viscElemsCoord, viscBlowing, reg.elemsCoord[:,:3])
-        else:
-            raise RuntimeError('Incorrect number of dimensions', self.cfg['nDim'])
+        # from matplotlib import pyplot as plt
+        # plt.figure()
+        # plt.plot(xv_up[:,0], blwv_up, label='visc up')
+        # plt.plot(xi_up[:,0], blwi_up, '--', label='inv up')
+        # plt.legend(frameon=False)
+        # plt.draw()
+        # plt.figure()
+        # plt.plot(xv_lw[:,0], blwv_lw, label='visc lw')
+        # plt.plot(xi_lw[:,0], blwi_lw, '--', label='inv lw')
+        # plt.legend(frameon=False)
+        # plt.draw()
+        # plt.figure()
+        # plt.plot(xv_wk[:,0], blwv_wk, label='visc wk')
+        # plt.plot(xi_wk[:,0], blwi_wk, '--', label='inv wk')
+        # plt.legend(frameon=False)
+        # plt.draw()
+        # plt.title('Blowing velocity')
+        # plt.show()
 
     def __rbfToSection(self, x, var, xs):
         if np.all(var == var[0]):
diff --git a/blast/tests/apiTest.py b/blast/tests/apiTest.py
index 9b38c095c1407fb3f507732d4de1a68020ce5c0d..7411c62703adfea86d3373ff6b0899787273b6dd 100644
--- a/blast/tests/apiTest.py
+++ b/blast/tests/apiTest.py
@@ -87,7 +87,7 @@ def cfgBlast(verb):
         'CFL0' : 1,                     # Inital CFL number of the calculation
         'Verb': verb,                   # Verbosity level of the solver
         'couplIter': 100,               # Maximum number of iterations
-        'couplTol' : 1e-2,              # Tolerance of the VII methodology
+        'couplTol' : 1e-4,              # Tolerance of the VII methodology
         'iterPrint': 20,                # int, number of iterations between outputs
         'resetInv' : True,              # bool, flag to reset the inviscid calculation at every iteration.
         'sections' : [0],               # List of sections for boundary layer calculation
diff --git a/blast/tests/dart/naca0012_rbf_2D.py b/blast/tests/dart/naca0012_rbf_2D.py
new file mode 100644
index 0000000000000000000000000000000000000000..29c53774322c1af4da25fa485afe17550f5b7bb4
--- /dev/null
+++ b/blast/tests/dart/naca0012_rbf_2D.py
@@ -0,0 +1,198 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2022 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.
+
+
+# @author Paul Dechamps
+# @date 2022
+# Test the blast implementation. The test case is a compressible attached transitional flow.
+# Tested functionalities;
+# - Time integration.
+# - Two time-marching techniques (agressive and safe).
+# - Transition routines.
+# - Forced transition.
+# - Compressible flow routines.
+# - Laminar and turbulent flow.
+#
+# Untested functionalities.
+# - Separation routines.
+# - Multiple failure case at one iteration.
+# - Response to unconverged inviscid solver.
+
+# Imports.
+
+import blast.utils as viscUtils
+import numpy as np
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def cfgInviscid(nthrds, verb):
+    import os.path
+    # Parameters
+    return {
+    # Options
+    'scenario' : 'aerodynamic',
+    'task' : 'analysis',
+    'Threads' : nthrds, # number of threads
+    'Verb' : verb, # verbosity
+    # Model (geometry or mesh)
+    'File' : os.path.dirname(os.path.abspath(__file__)) + '/../../models/dart/n0012.geo', # Input file containing the model
+    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 8.33, 'msTe' : 0.01, 'msLe' : 0.001}, # parameters for input file model
+    'Dim' : 2, # problem dimension
+    'Format' : 'gmsh', # save format (vtk or gmsh)
+    # Markers
+    'Fluid' : 'field', # name of physical group containing the fluid
+    'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+    'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
+    'Wake' : 'wake', # LIST of names of physical group containing the wake
+    'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
+    'Te' : 'te', # LIST of names of physical group containing the trailing edge
+    'dbc' : True,
+    'Upstream' : 'upstream',
+    # Freestream
+    'M_inf' : 0.2, # freestream Mach number
+    'AoA' : 2., # freestream angle of attack
+    # Geometry
+    'S_ref' : 1., # reference surface length
+    'c_ref' : 1., # reference chord length
+    'x_ref' : .25, # reference point for moment computation (x)
+    'y_ref' : 0.0, # reference point for moment computation (y)
+    'z_ref' : 0.0, # reference point for moment computation (z)
+    # Numerical
+    'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
+    'G_fill' : 2, # fill-in factor for GMRES preconditioner
+    'G_tol' : 1e-5, # tolerance for GMRES
+    'G_restart' : 50, # restart for GMRES
+    'Rel_tol' : 1e-6, # relative tolerance on solver residual
+    'Abs_tol' : 1e-8, # absolute tolerance on solver residual
+    'Max_it' : 20 # solver maximum number of iterations
+    }
+
+def cfgBlast(verb):
+    return {
+        'Re' : 1e7,                 # Freestream Reynolds number
+        'Minf' : 0.2,               # Freestream Mach number (used for the computation of the time step only)
+        'CFL0' : 1,                 # Inital CFL number of the calculation
+        'Verb': verb,               # Verbosity level of the solver
+        'couplIter': 200,            # Maximum number of iterations
+        'couplTol' : 1e-4,          # Tolerance of the VII methodology
+        'iterPrint': 1,             # int, number of iterations between outputs
+        'resetInv' : True,          # bool, flag to reset the inviscid calculation at every iteration.
+        'sections' : [0],           # List of sections for boundary layer calculation
+        'xtrF' : [None, 0.4],       # Forced transition location
+        'interpolator' : 'Rbf',      # Interpolator for the coupling
+        'neighbors': 10,
+        'rbftype': 'linear',
+        'smoothing': 1e-8,
+        'degree': 0.,
+        'span' : 1.,
+        'nPoints' : 170,
+    }
+
+def main():
+    # Timer.
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    args = parseargs()
+    icfg = cfgInviscid(args.k, args.verb)
+    vcfg = cfgBlast(args.verb)
+
+    tms['pre'].start()
+    coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
+    tms['pre'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    aeroCoeffs = coupler.run()
+    tms['solver'].stop()
+
+    # Display results.
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
+
+     # Write results to file.
+    vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all')
+    tms['total'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    print('SOLVERS statistics')
+    print(coupler.tms)
+    
+    # Test solution
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    tests.add(CTest('Cl', isol.getCl(), 0.230, 5e-2)) # XFOIL 0.2325
+    tests.add(CTest('Cd wake', vsol.Cdt, 0.0056, 1e-3, forceabs=True)) # XFOIL 0.00531
+    tests.add(CTest('Cd integral', vsol.Cdf + isol.getCd(), 0.0057, 1e-3, forceabs=True)) # XFOIL 0.00531
+    tests.add(CTest('Cdf', vsol.Cdf, 0.0046, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
+    tests.add(CTest('xtr_top', vsol.getAverageTransition(0), 0.195, 0.03, forceabs=True)) # XFOIL 0.1877
+    tests.add(CTest('xtr_bot', vsol.getAverageTransition(1), 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
+    tests.add(CTest('error dStar', np.linalg.norm(vSolution[0]['deltaStar'] - vSolution[0]['theta']*vSolution[0]['H']), 0., 0., forceabs=True))
+    tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 30, 0, forceabs=True))
+    tests.run()
+
+    # Show results
+    if not args.nogui:
+        iCp = viscUtils.read('Cp_inviscid.dat')
+        vCp = viscUtils.read('Cp_viscous.dat')
+        xfoilCp = viscUtils.read('../../blast/models/references/cpXfoil_n0012.dat', delim=None, skip = 1)
+        xfoilCpInv = viscUtils.read('../../blast/models/references/cpXfoilInv_n0012.dat', delim=None, skip = 1)
+        plotcp = {'curves': [np.column_stack((vCp[:,0], vCp[:,3])),
+                            np.column_stack((xfoilCp[:,0], xfoilCp[:,1])),
+                            np.column_stack((iCp[:,0], iCp[:,3])),
+                            np.column_stack((xfoilCpInv[:,0], xfoilCpInv[:,1]))],
+                    'labels': ['Blaster (VII)', 'XFOIL (VII)', 'DART (inviscid)', 'XFOIL (inviscid)'],
+                    'lw': [3, 3, 2, 2],
+                    'color': ['firebrick', 'darkblue', 'firebrick', 'darkblue'],
+                    'ls': ['-', '-', '--', '--'],
+                    'reverse': True,
+                    'xlim':[0, 1],
+                    'yreverse': True,
+                    'legend': True,
+                    'xlabel': '$x/c$',
+                    'ylabel': '$c_p$'
+                    }
+        viscUtils.plot(plotcp)
+
+        xfoilBl = viscUtils.read('../../blast/models/references/blXfoil_n0012.dat', delim=None, skip = 1)
+        plotcf = {'curves': [np.column_stack((vSolution[0]['x'], vSolution[0]['cf']*vSolution[0]['ue']*vSolution[0]['ue'])),
+                            np.column_stack((xfoilBl[:,1], xfoilBl[:,6]))],
+                'labels': ['Blaster', 'XFOIL'],
+                'lw': [3, 3],
+                'color': ['firebrick', 'darkblue'],
+                'ls': ['-', '-'],
+                'reverse': True,
+                'xlim':[0, 1],
+                'legend': True,
+                'xlabel': '$x/c$',
+                'ylabel': '$c_f$'
+                    }
+        viscUtils.plot(plotcf)
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()
diff --git a/blast/tests/dart/onera_3D.py b/blast/tests/dart/onera_3D.py
deleted file mode 100644
index 2ad42d166fd0c67413dcff7a633233836ee4d034..0000000000000000000000000000000000000000
--- a/blast/tests/dart/onera_3D.py
+++ /dev/null
@@ -1,157 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2022 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.
-
-
-# @author Paul Dechamps
-# @date 2022
-# Test the blast implementation on the 3D ONERA M6 wing.
-
-# Imports.
-
-import blast.utils as viscUtils
-import numpy as np
-import os.path
-
-from fwk.wutils import parseargs
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-import math
-
-def cfgInviscid(nthrds, verb):
-    # Parameters
-    return {
-    # Options
-    'scenario' : 'aerodynamic',
-    'task' : 'analysis',
-    'Threads' : nthrds, # number of threads
-    'Verb' : verb, # verbosity
-    # Model (geometry or mesh)
-    'File' : os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + '/models/dart/oneraM6.geo', # Input file containing the model
-    'Pars' : {'lgt' : 3.0, 'wdt' : 3.0, 'hgt' : 3.0, 'msLeRt' : 0.004, 'msTeRt' : 0.008, 'msLeTp' : 0.002, 'msTeTp' : 0.008, 'msF' : 0.6}, # parameters for input file model
-    'Dim' : 3, # problem dimension
-    'Format' : 'vtk', # save format (vtk or gmsh)
-    # Markers
-    'Fluid' : 'field', # name of physical group containing the fluid
-    'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
-    'Wings' : ['wing'], # LIST of names of physical groups containing the lifting surface boundary
-    'Wakes' : ['wake'], # LIST of names of physical group containing the wake
-    'WakeTips' : ['wakeTip'], # LIST of names of physical group containing the edge of the wake
-    'Tes' : ['te'], # LIST of names of physical group containing the trailing edge
-    'dbc' : True,
-    'Upstream' : 'upstream',
-    # Freestream
-    'M_inf' : 0.839,     # freestream Mach number
-    'AoA' : 3.06,        # freestream angle of attack
-    # Geometry
-    'S_ref' : 0.7528,    # reference surface length
-    'c_ref' : 0.64,      # reference chord length
-    'x_ref' : 0.0,       # reference point for moment computation (x)
-    'y_ref' : 0.0,       # reference point for moment computation (y)
-    'z_ref' : 0.0,       # reference point for moment computation (z)
-    # Numerical
-    'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
-    'G_fill' : 2,        # fill-in factor for GMRES preconditioner
-    'G_tol' : 1e-5,      # tolerance for GMRES
-    'G_restart' : 50,    # restart for GMRES
-    'Rel_tol' : 1e-6,    # relative tolerance on solver residual
-    'Abs_tol' : 1e-8,    # absolute tolerance on solver residual
-    'Max_it' : 50        # solver maximum number of iterations
-    }
-
-def cfgBlast(verb):
-    return {
-        'Re' : 11.72*10e6,  # Freestream Reynolds number
-        'Minf' : 0.839,     # Freestream Mach number (used for the computation of the time step only)
-        'CFL0' : 1,         # Inital CFL number of the calculation
-
-        'sections' : np.linspace(0.01, 1.1, 20),
-        'nPoints': 150,
-        'writeSections': [0.20, 0.44, 0.80],
-        'Sym':[0.],
-        'span':1.196,
-        'interpolator': 'Rbf',
-        'rbftype': 'linear',
-        'smoothing': 1e-8,
-        'degree': 0,
-        'neighbors': 10,
-        'saveTag': 5,
-
-        'Verb': verb,       # Verbosity level of the solver
-        'couplIter': 50,    # Maximum number of iterations
-        'couplTol' : 5e-4,  # Tolerance of the VII methodology
-        'iterPrint': 5,     # int, number of iterations between outputs
-        'resetInv' : True,  # bool, flag to reset the inviscid calculation at every iteration.
-        'xtrF' : [0.01, 0.01],# Forced transition location
-        'nDim' : 3
-    }
-
-def main():
-    # Timer.
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    args = parseargs()
-    icfg = cfgInviscid(args.k, args.verb)
-    vcfg = cfgBlast(args.verb)
-
-    parsViscous = {'nLe': 20, 'nTe': 8, 'nMid': 40, 'nSpan': 60, 'nWake': 20}
-    #vMsh = viscUtils.mesh(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/onera_visc.geo', parsViscous)
-    #vcfg['vMsh'] = vMsh
-
-    tms['pre'].start()
-    coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
-    tms['pre'].stop()
-
-    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-    tms['solver'].start()
-    aeroCoeffs = coupler.run()
-    tms['solver'].stop()
-
-    # Display results.
-    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
-    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
-    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
-
-     # Write results to file.
-    vSolution = viscUtils.getSolution(isol.sec, write=True, toW='all', sfx='onera')
-
-    # Save pressure coefficient
-    isol.save(sfx='_viscous')
-    tms['total'].stop()
-
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    print('SOLVERS statistics')
-    print(coupler.tms)
-    
-    # Test solution
-    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
-    tests = CTests()
-    tests.add(CTest('Cl', isol.getCl(), 0.279, 5e-2))
-    tests.add(CTest('Cd wake', vsol.Cdt, 0.00471, 1e-3, forceabs=True))
-    tests.add(CTest('Cd int', isol.getCd() + vsol.Cdf, 0.01519, 1e-3, forceabs=True))
-    tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 13, 0, forceabs=True))
-    tests.run()
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/blast/tests/dart/rae2822_3D.py b/blast/tests/dart/rae2822_3D.py
index 61167b6dee02949d67f2eca4cd9cb8ae81829761..4c6de7dff8cd56205253036976b6f3a344c805d4 100644
--- a/blast/tests/dart/rae2822_3D.py
+++ b/blast/tests/dart/rae2822_3D.py
@@ -55,7 +55,7 @@ def cfgInviscid(nthrds, verb):
     'Verb' : verb, # verbosity
     # Model (geometry or mesh)
     'File' : os.path.abspath(os.path.join(os.path.abspath(__file__), '../../../models/dart/rae_3.geo')), # Input file containing the model
-    'Pars' : {'spn' : 1.0, 'lgt' : 6.0, 'wdt' : 3.0, 'hgt' : 6.0, 'msN' : 0.02, 'msF' : 1}, # parameters for input file model
+    'Pars' : {'spn' : 1.0, 'lgt' : 6.0, 'wdt' : 3.0, 'hgt' : 6.0, 'msN' : 0.03, 'msF' : 1.2}, # parameters for input file model
     'Dim' : 3, # problem dimension
     'Format' : 'vtk', # save format (vtk or gmsh)
     # Markers
@@ -93,7 +93,8 @@ def cfgBlast(verb):
         'Minf' : 0.8,     # Freestream Mach number (used for the computation of the time step only)
         'CFL0' : 1,       # Inital CFL number of the calculation
 
-        'sections' : np.linspace(0.01, 0.95, 3),
+        'sections' : np.linspace(0.01, 0.95, 50),
+        'nPoints' : 500,
         'writeSections': [0.2, 0.4, 0.6, 0.8, 1.0],
         'Sym':[0.],
         'span': 1.,
@@ -109,7 +110,7 @@ def cfgBlast(verb):
         'couplTol' : 5e-2,  # Tolerance of the VII methodology
         'iterPrint': 1,     # int, number of iterations between outputs
         'resetInv' : False, # bool, flag to reset the inviscid calculation at every iteration.
-        'xtrF' : [0., 0.],  # Forced transition location
+        'xtrF' : [0.01, 0.01],  # Forced transition location
         'nDim' : 3
     }
 
@@ -122,21 +123,13 @@ def main():
     icfg = cfgInviscid(args.k, args.verb)
     vcfg = cfgBlast(args.verb)
 
-    parsViscous = {'spn': icfg['Pars']['spn'], 'lgt': icfg['Pars']['lgt'],
-                   'nLe': 25, 'nMid': 50, 'nTe': 10, 'nSpan': 60, 'nWake': 25,
-                   'progLe': 1.1, 'progMid': 1.0,  'progTe': 1.0, 'progSpan': 1.0, 'progWake': 1.15}
-    
-    vMeshFile = os.path.abspath(os.path.join(os.path.abspath(__file__), '../../../models/dart/rae_3_visc.geo'))
-    vMsh = viscUtils.mesh(vMeshFile, parsViscous)
-    vcfg['vMsh'] = vMsh
-
     tms['pre'].start()
     coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
     tms['pre'].stop()
 
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
     tms['solver'].start()
-    aeroCoeffs = coupler.run(write=False)
+    aeroCoeffs = coupler.run(write=True)
     tms['solver'].stop()
 
     # Display results.
@@ -160,10 +153,10 @@ def main():
     # Test solution
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
-    tests.add(CTest('Cl', isol.getCl(), 0.135, 5e-2))
-    tests.add(CTest('Cd wake', vsol.Cdt, 0.0074, 1e-3, forceabs=True))
-    tests.add(CTest('Cd integral', vsol.Cdf + isol.getCd(), 0.0140, 1e-3, forceabs=True))
-    tests.add(CTest('Cdf', vsol.Cdf, 0.0061, 1e-3, forceabs=True))
+    tests.add(CTest('Cl', isol.getCl(), 0.10, 5e-2))
+    tests.add(CTest('Cd wake', vsol.Cdt, 0.0044, 1e-3, forceabs=True))
+    tests.add(CTest('Cd integral', vsol.Cdf + isol.getCd(), 0.0128, 1e-3, forceabs=True))
+    tests.add(CTest('Cdf', vsol.Cdf, 0.0063, 1e-3, forceabs=True))
     tests.add(CTest('Iterations', len(aeroCoeffs), 3, 0, forceabs=True))
     tests.run()