diff --git a/blast/api/core.py b/blast/api/core.py
index 4d6b826bc9362cdd1d27e7e9f145f51be1e6601b..5811adfd4c3816d07972fa56f887b08839c7b4b0 100644
--- a/blast/api/core.py
+++ b/blast/api/core.py
@@ -57,8 +57,6 @@ def initBlast(cfg, icfg, iSolverName='DART'):
 
     if iSolverName == 'DART':
         from blast.interfaces.dart.DartInterface import DartInterface as interface
-        from dart.api.core import initDart
-        iSolverObjects = initDart(icfg, scenario=icfg['scenario'], task=icfg['task'], viscous = 1)
 
     # Check viscous solver parameters.
     if 'Re' in cfg and cfg['Re'] > 0:
@@ -105,13 +103,12 @@ def initBlast(cfg, icfg, iSolverName='DART'):
         __resetInv = False
 
     # Viscous solver object.
-    vSolver = blast.Driver(_Re, _Minf, _CFL0, _nSections, _xtrF[0], _xtrF[1], _span, verbose =_verbose)
+    vSolver = blast.Driver(_Re, _Minf, _CFL0, _nSections, _xtrF[0], _xtrF[1], _span, _verbose =_verbose)
     # Solvers interface.
-    solversAPI = interface(iSolverObjects['sol'], vSolver, [iSolverObjects['blwb'], iSolverObjects['blww']])
+    solversAPI = interface(icfg, vSolver, cfg)
     # Coupler
     coupler = Coupler(solversAPI, vSolver, _maxCouplIter = __couplIter, _couplTol = __couplTol, _iterPrint = __iterPrint, _resetInv = __resetInv)
 
     return {'coupler': coupler,
             'solversAPI': solversAPI,
-            'vSolver': vSolver,
-            'iSolverObjects': iSolverObjects}
+            'vSolver': vSolver}
diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index cdfaa70c6f281d14020ee92d23be253997433d1b..aab4547a575765ff01cce87269f93a4a7c2087d4 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/blast/interfaces/DSolversInterface.py
@@ -22,7 +22,7 @@ class SolversInterface:
         elif self.cfg['interpolator'] == 'Rbf':
             from blast.interfaces.interpolators.DRbfInterpolator import RbfInterpolator as interp
             print('RBF interpolator initialized.')
-            print('Loading viscous mesh... ')
+            print('Loading viscous mesh... ', end='')
             if self.cfg['nDim'] == 2:
                 self.getVAirfoil()
             elif self.cfg['nDim'] == 3:
@@ -32,7 +32,7 @@ class SolversInterface:
             raise RuntimeError('Incorrect interpolator specified in config.')
         self.interpolator = interp(self.cfg)
 
-        print('Loading inviscid mesh... ')
+        print('Loading inviscid mesh... ', end='')
         self.setMesh()
         print('done.\n')
 
@@ -50,42 +50,6 @@ class SolversInterface:
         """
         self.interpolator.inviscidToViscous(self.iBnd, self.vBnd, couplIter)
         for iSec in range(self.cfg['nSections']):
-            """print('')
-            print('// Section ', iSec)
-            for i in range(len(self.vBnd[iSec][0].nodesCoord[:,0])):
-                print('Point({0:>5.0f}) = ({1:>8.6f}, {2:>8.6f}, {3:>8.6f});'.format(i+(iSec+7)*10000, self.vBnd[iSec][0].nodesCoord[i,0], self.vBnd[iSec][0].nodesCoord[i,2], self.vBnd[iSec][0].nodesCoord[i,1]))"""
-            """from matplotlib import pyplot as plt
-            plt.plot(self.vBnd[iSec][0].getXUpper(), self.vBnd[iSec][0].getYUpper(), 'o-')
-            plt.plot(self.vBnd[iSec][0].getXLower(), self.vBnd[iSec][0].getYLower(), 'x-')
-            plt.plot(self.vBnd[iSec][1].nodesCoord[:,0], self.vBnd[iSec][1].nodesCoord[:,1], 'x-')
-            plt.show()
-            plt.plot(self.vBnd[0][0].getXUpper(), self.vBnd[0][0].getVelocityXUpper())
-            plt.plot(self.vBnd[0][0].getXLower(), self.vBnd[0][0].getVelocityXLower())
-            #plt.plot(self.vBnd[0][1].nodesCoord[:,0], self.vBnd[0][1].V[:,0])
-            plt.title('vx')
-            plt.show()
-            plt.plot(self.vBnd[0][0].getXUpper(), self.vBnd[0][0].getVelocityYUpper())
-            plt.plot(self.vBnd[0][0].getXLower(), self.vBnd[0][0].getVelocityYLower())
-            plt.plot(self.vBnd[0][1].nodesCoord[:,0], self.vBnd[0][1].V[:,1])
-            plt.title('vy')
-            plt.show()
-            plt.plot(self.vBnd[0][0].getXUpper(), self.vBnd[0][0].getVelocityZUpper())
-            plt.plot(self.vBnd[0][0].getXLower(), self.vBnd[0][0].getVelocityZLower())
-            plt.plot(self.vBnd[0][1].nodesCoord[:,0], self.vBnd[0][1].V[:,2])
-            plt.title('vz')
-            plt.show()
-            plt.plot(self.vBnd[0][0].getXUpper(), self.vBnd[0][0].getMachUpper())
-            plt.plot(self.vBnd[0][0].getXLower(), self.vBnd[0][0].getMachLower())
-            plt.plot(self.vBnd[0][1].nodesCoord[:,0], self.vBnd[0][1].M)
-            plt.title('Me')
-            plt.show()
-            plt.plot(self.vBnd[0][0].getXUpper(), self.vBnd[0][0].getRhoUpper())
-            plt.plot(self.vBnd[0][0].getXLower(), self.vBnd[0][0].getRhoLower())
-            plt.plot(self.vBnd[0][1].nodesCoord[:,0], self.vBnd[0][1].Rho)
-            plt.title('Rho')
-            plt.show()
-            quit()"""
-
             ## Mesh
             if (self.vBnd[iSec][0].isMeshChanged()):
                 for reg in self.vBnd[iSec]:
@@ -218,7 +182,7 @@ class SolversInterface:
         self.cfg['EffSections'] = np.empty(0)
 
         for iReg, reg in enumerate(data):
-            for iy, y in enumerate(self.cfg['Sections']):
+            for iy, y in enumerate(self.cfg['sections']):
                 # Find nearest value in the mesh
                 try:
                     idx = (np.abs(reg[:,2]-y).argmin())
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index 321fec380bfa925a2e31289f21ee867d9ad906c2..10c0f11e3ce2d6c5cd2144f3c9452b875d860f35 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/blast/interfaces/dart/DartInterface.py
@@ -36,6 +36,9 @@ class DartInterface(SolversInterface):
         except:
             print(ccolors.ANSI_RED, 'Could not load DART. Make sure it is installed.', ccolors.ANSI_RESET)
             raise RuntimeError('Import error')
+
+        if 'iMsh' not in _cfg:
+            _cfg['iMsh'] = self.blw
         SolversInterface.__init__(self, vSolver, _cfg)
     
     def run(self):
@@ -298,8 +301,6 @@ class DartInterface(SolversInterface):
             raise RuntimeError('Could not identify the lower side in the inviscid mesh.')
         """
 
-        print('ddd')
-
         for iRegion in range(len(self.cfg['iMsh'])):
             for e in self.cfg['iMsh'][iRegion].tag.elems:
                 # Compute cg position
@@ -315,12 +316,7 @@ class DartInterface(SolversInterface):
                 cg_pt = np.hstack((cg_pt, e.no))
                 cg[iRegion] = np.vstack((cg[iRegion], cg_pt))
 
-        print('cccccc')
         self.ilwIdx = []
-        """for idx in range(len(data[0][:,3])):
-            if data[0][idx, 3] in lwRows:
-                self.ilwIdx.append(idx)"""
-
         if 'Sym' in self.cfg:
             for s in self.cfg['Sym']:
                 for iReg in range(len(data)):
diff --git a/blast/interfaces/interpolators/DMatchingInterpolator.py b/blast/interfaces/interpolators/DMatchingInterpolator.py
index e286024b6f3eec8c006f1efa6ffde482f18e02f8..a8018902977847edc7e22b7cec598d2cba4772e5 100644
--- a/blast/interfaces/interpolators/DMatchingInterpolator.py
+++ b/blast/interfaces/interpolators/DMatchingInterpolator.py
@@ -10,14 +10,11 @@ class MatchingInterpolator:
             for iReg in range(len(iDict)):
                 vDict[0][iReg].updateVars(iDict[iReg].nodesCoord, iDict[iReg].V, iDict[iReg].M, iDict[iReg].Rho)
         elif self.cfg['nDim'] == 3:
-            for iSec, ysec in enumerate(self.cfg['Sections']):
+            for iSec, ysec in enumerate(self.cfg['sections']):
                 for iReg in range(2):
-                    print(iSec)
                     print(iDict[iReg].nodesCoord[iDict[iReg].nodesCoord[:,1] == ysec])
                     print(iDict[iReg].V[iDict[iReg].nodesCoord[:,1] == ysec])
                     vDict[iSec][iReg].updateVars(iDict[iReg].nodesCoord[iDict[iReg].nodesCoord[:,1] == ysec], iDict[iReg].V[iDict[iReg].nodesCoord[:,1] == ysec], iDict[iReg].Rho[iDict[iReg].nodesCoord[:,1] == ysec])
-                    
-                    
 
     def viscousToInviscid(self, iDict, vDict):
         if self.cfg['nDim'] == 2:
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index ab58b07e8236e78fe81d6d26b0a07e7bf54aabb2..08a2e328651df0cee65513e5d7fc91a7a3061421 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -528,6 +528,8 @@ int Driver::outputStatus(double maxMach) const
  */
 void Driver::setGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z)
 {
+    if (iSec > sections.size() || iRegion > sections[iSec].size())
+        throw std::runtime_error("blast::Driver Wrong section or region id.");
     sections[iSec][iRegion]->setMesh(_x, _y, _z);
 }
 
diff --git a/blast/tests/dart/blidart.py b/blast/tests/dart/blidart.py
index a1549e1cfca57e8b32da4c7a2da9e91beb7fc7a2..16b164d71c6d4ea36bf747d0c273d972b8b98b21 100644
--- a/blast/tests/dart/blidart.py
+++ b/blast/tests/dart/blidart.py
@@ -148,7 +148,7 @@ def main():
     tests.add(CTest('Cdf', vSolution['Cdf'], 0.00465, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
     tests.add(CTest('xtr_top', vSolution['xtrT'], 0.20, 0.03, forceabs=True)) # XFOIL 0.1877
     tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
-    tests.add(CTest('Iterations', len(aeroCoeffs), 18, 0, forceabs=True))
+    tests.add(CTest('Iterations', len(aeroCoeffs), 17, 0, forceabs=True))
     tests.run()
 
     # Show results
diff --git a/blast/utils.py b/blast/utils.py
index b84f57609df78fe8eb0b24dd956e7f9453ee179c..b07f83c3a79fa1ef1d910b132ea96cf16ea407f1 100644
--- a/blast/utils.py
+++ b/blast/utils.py
@@ -64,8 +64,11 @@ def initBlast(iconfig, vconfig, iSolver='DART'):
     - iSolverAPI: api interface of the inviscid solver selected with 'iSolver'.
     - vSolver: blast::Driver class.
     """
+
+    if 'nSections' not in vconfig:
+        vconfig['nSections'] = len(vconfig['sections'])
     # Viscous solver
-    vSolver = initBL(vconfig['Re'], vconfig['Minf'], vconfig['CFL0'], 1, xtrF=vconfig['xtrF'])
+    vSolver = initBL(vconfig['Re'], vconfig['Minf'], vconfig['CFL0'], vconfig['nSections'], xtrF=vconfig['xtrF'])
 
     # Inviscid solver
     if iSolver == 'DART':