diff --git a/blast/interfaces/blSolversInterface.py b/blast/interfaces/blSolversInterface.py
index 5f3c4cf849758ae387a016b6fc32bd796328846c..ced73ad27bf6d8725957a3d16a415d69b76852d7 100644
--- a/blast/interfaces/blSolversInterface.py
+++ b/blast/interfaces/blSolversInterface.py
@@ -22,7 +22,9 @@ import numpy as np
 import blast
 import tbox
 
-from blast.interfaces.blDataStructure import Group
+# from blast.interfaces.blDataStructure import Group
+from blast.interfaces.blDataInviscid import InviscidData
+from blast.interfaces.blDataViscous import ViscousData
 from scipy.interpolate import griddata
 from scipy.interpolate import interp1d
 
@@ -90,17 +92,21 @@ class SolversInterface:
             print(f'-Creating sections using {self.cfg["genSections"]}', end='...')
             self._getSections()
             print('done')
-
-        # External quantities
+        
         self.deltaStarExt = None
-        self.xxExt        = None
-        self.vtExt        = None
+        self.xxExt = None
+        self.vtExt = None
 
-        self.vSolver = None
-        # Flags
         self.vinit = False
+        # Flags
+        self.vSolver = None
         self.hasvsol = False
 
+        # Check for restart file
+        if "restart_solution" in self.cfg and self.cfg["restart_solution"] == True:
+            print('Loading restart file', end='...')
+            self._load_restart(self.cfg["restart_inputName"], self.vBnd)
+            print('done')
         self._printSetup()
 
     ########################
@@ -125,33 +131,41 @@ class SolversInterface:
             couplIter: int
                 Current iteration of the coupling loop.
         """
-        if couplIter != 0:
-            for ibody, body in enumerate(self.vSolver.bodies):
-                for isec, section in enumerate(body.sections):
-                    for iReg in range(len(self.vBnd[ibody][isec])):
-                        pystruct = self.vBnd[ibody][isec][iReg]
-
-                        if 'vWake' in pystruct.name or pystruct.type == 'fuselage':
-                            wakeId = section.getRegionIdByName("wake")
-                            if pystruct.blowingVel.shape[0] != len(section.regions[wakeId].blowingVelocity):
-                                raise RuntimeError('Blowing velocity size mismatch')
-                            pystruct.blowingVel = np.asarray(section.regions[wakeId].blowingVelocity)
-                        elif 'vAirfoil' in pystruct.name:
-                            upperId = section.getRegionIdByName("upper")
-                            lowerId = section.getRegionIdByName("lower")
-                            if pystruct.blowingVel.shape[0] != len(np.concatenate((np.flip(np.asarray(section.regions[upperId].blowingVelocity)),
-                                                                            np.asarray(section.regions[lowerId].blowingVelocity)))):
-                                raise RuntimeError('Blowing velocity size mismatch')
-                            pystruct.blowingVel = np.concatenate((np.flip(np.asarray(section.regions[upperId].blowingVelocity)),
-                                                            np.asarray(section.regions[lowerId].blowingVelocity)))
+        if couplIter == 0:
+            return
+        for ibody, body in enumerate(self.vSolver.bodies):
+            for isec, section in enumerate(body.sections):
+                for iReg, pystruct in enumerate(self.vBnd[ibody][isec]):
 
-                    for i, reg in enumerate(section.regions):
-                        loc = np.zeros(reg.getnNodes())
-                        for inod in range(reg.getnNodes()):
-                            loc[inod] = reg.nodes[inod].xi
-                        self.xxExt[ibody][isec][i] = np.asarray(loc)
-                        self.deltaStarExt[ibody][isec][i] = np.asarray(reg.deltaStar)
-                        self.vtExt[ibody][isec][i] = np.asarray(reg.vt)
+                    name, ptype = pystruct.getName(), pystruct.getType()
+                    blowingVelocity = pystruct.getBlowingVelocity()
+
+                    if 'vWake' in name or ptype == 'fuselage':
+                        wakeId = section.getRegionIdByName("wake")
+                        wakeVel = np.asarray(section.regions[wakeId].blowingVelocity)
+
+                        if blowingVelocity.shape[0] != wakeVel.shape[0]:
+                            raise RuntimeError('Blowing velocity size mismatch')
+
+                        pystruct.setBlowingVelocity(wakeVel)
+
+                    elif 'vAirfoil' in name:
+                        upperId, lowerId = section.getRegionIdByName("upper"), section.getRegionIdByName("lower")
+                        upperVel = np.asarray(section.regions[upperId].blowingVelocity)
+                        lowerVel = np.asarray(section.regions[lowerId].blowingVelocity)
+
+                        combinedVel = np.concatenate((upperVel[::-1], lowerVel))
+
+                        if blowingVelocity.shape[0] != combinedVel.shape[0]:
+                            raise RuntimeError('Blowing velocity size mismatch')
+
+                        pystruct.setBlowingVelocity(combinedVel)
+
+                for i, reg in enumerate(section.regions):
+                    nodes_xi = np.array([node.xi for node in reg.nodes], dtype=np.float64)
+                    self.xxExt[ibody][isec][i] = nodes_xi
+                    self.deltaStarExt[ibody][isec][i] = np.array(reg.deltaStar, dtype=np.float64)
+                    self.vtExt[ibody][isec][i] = np.array(reg.vt, dtype=np.float64)
 
     def interpolate(self, dir:str):
         """Interpolate solution between inviscid and viscous meshes.
@@ -162,34 +176,46 @@ class SolversInterface:
                 'i2v': inviscid to viscous.
                 'v2i': viscous to inviscid.
         """
-        if dir == 'i2v':
-            self.interpolator.inviscidToViscous(self.iBnd, self.vBnd)
-        elif dir == 'v2i':
-            self.interpolator.viscousToInviscid(self.iBnd, self.vBnd)
-        else:
-            raise RuntimeError('Invalid direction', dir)
+        interpolation_methods = {
+        "i2v": self.interpolator.inviscidToViscous,
+        "v2i": self.interpolator.viscousToInviscid
+        }
+        try:
+            interpolation_methods[dir](self.iBnd, self.vBnd)
+        except KeyError:
+            raise RuntimeError(f"Invalid direction: {dir}")
 
     def setViscousSolver(self):
         """Sets the actual conditions in vBnd to the viscous solver.
         """
+        region_map = {
+        "wake": 1,
+        "upper": 0,
+        "lower": 0,
+        "fuselage": 0
+        }
+
         for ibody, body in enumerate(self.vSolver.bodies):
             for isec, section in enumerate(body.sections):
                 for i, reg in enumerate(section.regions):
-                    if reg.getName() in ['wake']:
-                        iReg = 1
-                    elif reg.getName() in ['upper', 'lower', 'fuselage']:
-                        iReg = 0
-                    else:
-                        raise RuntimeError('Invalid region', reg.getName())
+                    reg_name = reg.getName()
+                    
+                    if reg_name not in region_map:
+                        raise RuntimeError(f"Invalid region: {reg_name}")
+
+                    iReg = region_map[reg_name]
                     pystruct = self.vBnd[ibody][isec][iReg]
-                    Mach = pystruct.getMach(reg.getName())
-                    Rho  = pystruct.getRho(reg.getName())
-                    Vt   = pystruct.getVt(reg.getName())
-                    for inod in range(reg.getnNodes()):
-                        reg.Me[inod] = Mach[inod]
-                        reg.rhoe[inod] = Rho[inod]
-                        reg.vt[inod] = Vt[inod]
 
+                    Mach, Rho, Vt = (
+                        pystruct.getMach(reg_name),
+                        pystruct.getDensity(reg_name),
+                        pystruct.getVt(reg_name),
+                    )
+
+                    for inod, (mach, rho, vt) in enumerate(zip(Mach, Rho, Vt)):
+                        reg.Me[inod] = mach
+                        reg.rhoe[inod] = rho
+                        reg.vt[inod] = vt
                         reg.deltaStarExt[inod] = self.deltaStarExt[ibody][isec][i][inod]
                         reg.nodes[inod].xiExt = self.xxExt[ibody][isec][i][inod]
                         reg.vtExt[inod] = self.vtExt[ibody][isec][i][inod]
@@ -210,9 +236,8 @@ class SolversInterface:
             if len(body.sections) != len(self.vBnd[ibody]):
                 raise RuntimeError(f'Expected {len(self.vBnd[ibody])} sections in viscous solver, got {len(body.sections)}')
             for isec, section in enumerate(body.sections):
-
-                # Check that if self.vBnd[ibody][isec][0].type is 'wing', section mus have upper and lower
-                if self.vBnd[ibody][isec][0].type == 'wing':
+                # Check that if self.vBnd[ibody][isec][0].getType() is 'wing', section mus have upper and lower
+                if self.vBnd[ibody][isec][0].getType() == 'wing':
                     hasupper = False
                     haslower = False
                     for reg in section.regions:
@@ -221,10 +246,10 @@ class SolversInterface:
                         elif reg.getName() == "lower":
                             haslower = True
                     if not hasupper or not haslower:
-                        raise RuntimeError(f'Expected upper and lower sections for type {self.vBnd[ibody][isec][0].type}')
+                        raise RuntimeError(f'Expected upper and lower sections for type {self.vBnd[ibody][isec][0].getType()}')
 
                 # Check that the fuselage has no wake
-                elif self.vBnd[ibody][isec][0].type == 'fuselage':
+                elif self.vBnd[ibody][isec][0].getType() == 'fuselage':
                     for reg in section.regions:
                         if reg.getName() == "wake":
                             raise RuntimeError('Fuselage should not have wake')
@@ -237,7 +262,7 @@ class SolversInterface:
                 bodyhaswake = body.hasWake()
                 pystructhaswake = False
                 for iReg, reg in enumerate(self.vBnd[ibody][isec]):
-                    if reg.type == 'wake':
+                    if reg.getType() == 'wake':
                         pystructhaswake = True
                         break
                 if bodyhaswake != pystructhaswake:
@@ -255,12 +280,36 @@ class SolversInterface:
         # Set mesh
         self._setViscousMesh()
 
-        self.xxExt = [[[np.zeros(0) for iReg in range(len(self.vSolver.bodies[ibody].sections[isec].regions))]\
-                for isec in range(len(self.vSolver.bodies[ibody].sections))] for ibody in range(self.getnBodies())]
-        self.deltaStarExt = [[[np.zeros(0) for iReg in range(len(self.vSolver.bodies[ibody].sections[isec].regions))]\
-                for isec in range(len(self.vBnd[ibody]))] for ibody in range(self.getnBodies())]
-        self.vtExt = [[[np.zeros(0) for iReg in range(len(self.vSolver.bodies[ibody].sections[isec].regions))]\
-                for isec in range(len(self.vSolver.bodies[ibody].sections))] for ibody in range(self.getnBodies())]
+        # Initialize interaction law quantities
+        if not self.vinit:
+            self._initializeInteractionLaw()
+        print('Viscous solver initialized')
+        return 0
+    
+    def _initializeInteractionLaw(self):
+        self.deltaStarExt = []
+        self.xxExt = []
+        self.vtExt = []
+        for ibody, body in enumerate(self.vBnd):
+            self.deltaStarExt.append([])
+            self.xxExt.append([])
+            self.vtExt.append([])
+            for isec, sec in enumerate(body):
+                self.deltaStarExt[ibody].append([])
+                self.xxExt[ibody].append([])
+                self.vtExt[ibody].append([])
+                for i, vreg in enumerate(sec):
+                    if vreg.getType() == 'fuselage' or vreg.getType() == 'wake':
+                        self.deltaStarExt[ibody][isec].append([]) # only one side
+                        self.xxExt[ibody][isec].append([])
+                        self.vtExt[ibody][isec].append([])
+                    elif vreg.getType() == 'wing':
+                        self.deltaStarExt[ibody][isec].append([]) # upper and lower
+                        self.deltaStarExt[ibody][isec].append([])
+                        self.xxExt[ibody][isec].append([])
+                        self.xxExt[ibody][isec].append([])
+                        self.vtExt[ibody][isec].append([])
+                        self.vtExt[ibody][isec].append([])
 
         # External variables
         if not self.vinit:
@@ -279,9 +328,7 @@ class SolversInterface:
                         self.xxExt[ibody][isec][i] = loc
                         self.deltaStarExt[ibody][isec][i] = np.zeros(reg.getnNodes())
                         self.vtExt[ibody][isec][i] = np.zeros(reg.getnNodes())
-        self.vinit = True
-        print('Viscous solver initialized')
-        return 0
+            self.vinit = True
 
     def _setViscousMesh(self):
         """Set the mesh from vBnd to the viscous solver.
@@ -313,16 +360,16 @@ class SolversInterface:
                     pystruct = self.vBnd[ibody][isec][iReg]
                     # Create node
                     nodes = blast.std_vector_blNode()
-                    for inod, nod in enumerate(pystruct.getNodesCoord(reg.getName())):
+                    for inod, nod in enumerate(pystruct.getNodes(reg.getName())):
                         pos = tbox.Vector3d(nod[xIdx], nod[yIdx], nod[zIdx])
-                        row = int(pystruct.getNodesRow(reg.getName())[inod]) if self.cfg['interpolator'] == 'Matching' else inod
+                        row = int(pystruct.getNodeConnectList(reg.getName())[inod]) if self.cfg['interpolator'] == 'Matching' else inod
                         nodes.push_back(blast.blNode(inod, pos, row))
                     # Set nodes
                     reg.setMesh(nodes)
                     # Set elements connectivity
                     connectElems = tbox.std_vector_int()
                     if self.cfg['interpolator'] == 'Matching':
-                        for val in pystruct.getConnectList(reg.getName()):
+                        for val in pystruct.getElmsConnectList(reg.getName()):
                             connectElems.push_back(int(val))
                     else:
                         for i in range(reg.getnElms()):
@@ -359,28 +406,23 @@ class SolversInterface:
         # Update the mesh in data structures
         for ibody in range(self.getnBodies()):
             for ireg, reg in enumerate(self.iBnd[ibody]):
-                for i, row in enumerate(reg.connectListNodes):
+                nodes = np.zeros((meshList[ibody][ireg].shape[0], 3), dtype=float)
+                elems = np.zeros((meshList[ibody][ireg].shape[0]-1, 3), dtype=float)
+                for i, row in enumerate(reg.getConnectListNodes()):
                     for idim in range(self.getnDim()):
-                        reg.nodesCoord[i, idim] = meshList[ibody][ireg][row, idim]
-
-        # Recompute elements data
-        for ibody in range(self.getnBodies()):
-            for ireg, reg in enumerate(self.iBnd[ibody]):
-                for i in range(len(reg.nodesCoord[:,0])-1):
+                        nodes[i, idim] = meshList[ibody][ireg][row, idim]
+                for i in range(reg.getnElms()-1):
                     for idim in range(self.getnDim()):
-                        reg.elemsCoord[i, idim] = (reg.nodesCoord[i, idim] + reg.nodesCoord[i+1, idim])/2
+                        elems[i, idim] = (nodes[i, idim] + nodes[i+1, idim])/2
+                reg.updateMesh(nodes, elems)
 
         if self.cfg['interpolator'] == 'Matching':
             for ibody in range(self.getnBodies()):
                 for iSec in range(len(self.vBnd)):
                     for iReg, vreg in enumerate(self.vBnd[ibody][iSec]):
                         ireg = self.iBnd[ibody][iReg]
-                        for i, row in enumerate(ireg.connectListNodes):
-                            for idim in range(self.getnDim()):
-                                vreg.nodesCoord[i, idim] = ireg.nodesCoord[i, idim]
-                        for i, row in enumerate(ireg.connectListElems):
-                            for idim in range(self.getnDim()):
-                                vreg.elemsCoord[i, idim] = ireg.elemsCoord[i, idim]
+                        vreg.updateMesh(ireg.getNodes('all'), ireg.getElms('all'))
+
         elif self.cfg['interpolator'] == 'Rbf':
             upNodes, lwNodes, upElms, lwElms = self.getUpperLowerSides()
             for ibody in range(self.getnBodies()):
@@ -479,8 +521,8 @@ class SolversInterface:
             if not self.iBnd[ibody][0].sidesIdentified:
                 raise RuntimeError(f'Cannot generate sections in auto for body {bodyName}: sides not identified.')
 
-            upper_nodes = self.iBnd[ibody][0].getUpperNodes().copy()
-            lower_nodes = self.iBnd[ibody][0].getLowerNodes().copy()
+            upper_nodes = self.iBnd[ibody][0].getNodes('upper_fromLE')
+            lower_nodes = self.iBnd[ibody][0].getNodes('lower_fromLE')
 
             # Remove points higher than span
             if ndim == 3:
@@ -557,8 +599,9 @@ class SolversInterface:
 
             # Wake
             wake = np.zeros((0,3))
-            for i in range(self.iBnd[ibody][1].nPoints):
-                wake = np.row_stack((wake, [self.iBnd[ibody][1].nodesCoord[i,0], self.iBnd[ibody][1].nodesCoord[i,1], self.iBnd[ibody][1].nodesCoord[i,2]]))
+            wakeNodes = self.iBnd[ibody][1].getNodes('all')
+            for i in range(self.iBnd[ibody][1].getnNodes()):
+                wake = np.row_stack((wake, [wakeNodes[i,0], wakeNodes[i,1], wakeNodes[i,2]]))
 
             interpolated_wake = []
 
@@ -602,12 +645,9 @@ class SolversInterface:
                     elems_coords_wake[isec][i, :] = (sec[i, :] + sec[i+1, :]) / 2
 
             for iy in range(len(sections_final)):
-                self.vBnd[ibody][iy][0].initStructures(sections_final[iy].shape[0], elems_coords_wing[iy].shape[0])
-                self.vBnd[ibody][iy][0].nodesCoord = sections_final[iy]
-                self.vBnd[ibody][iy][0].elemsCoord = elems_coords_wing[iy]
-                self.vBnd[ibody][iy][1].initStructures(wake_final[iy].shape[0], elems_coords_wake[iy].shape[0])
-                self.vBnd[ibody][iy][1].nodesCoord = wake_final[iy]
-                self.vBnd[ibody][iy][1].elemsCoord = elems_coords_wake[iy]
+                self.vBnd[ibody][iy][0].setMesh(sections_final[iy], elems_coords_wing[iy])
+                if len(self.vBnd[ibody][iy]) > 1:
+                    self.vBnd[ibody][iy][1].setMesh(wake_final[iy], elems_coords_wake[iy])
 
     def _getSectionsFromVTK(self):
 
@@ -669,15 +709,12 @@ class SolversInterface:
                 elems_coords_wing = np.empty((section.shape[0]-1, section.shape[1]), dtype=float, order='C')
                 for ielm in range(section.shape[0]-1):
                     elems_coords_wing[ielm, :] = (section[ielm, :] + section[ielm+1, :]) / 2
-
-                self.vBnd[ibody][isec][0].initStructures(section.shape[0], elems_coords_wing.shape[0])
-                self.vBnd[ibody][isec][0].nodesCoord = section
-                self.vBnd[ibody][isec][0].elemsCoord = elems_coords_wing
+                self.vBnd[ibody][isec][0].setMesh(section, elems_coords_wing)
 
                 # Wake
                 hasWake = False
                 for ibnd in self.iBnd[ibody]:
-                    if ibnd.type == 'wake':
+                    if ibnd.getType() == 'wake':
                         hasWake = True
                         break
 
@@ -696,9 +733,7 @@ class SolversInterface:
                     for ielm in range(wake_final.shape[0]-1):
                         elems_coords_wake[ielm, :] = (wake_final[ielm, :] + wake_final[ielm+1, :]) / 2
 
-                    self.vBnd[ibody][isec][1].initStructures(wake_final.shape[0], elems_coords_wake.shape[0])
-                    self.vBnd[ibody][isec][1].nodesCoord = wake_final
-                    self.vBnd[ibody][isec][1].elemsCoord = elems_coords_wake
+                    self.vBnd[ibody][isec][1].setMesh(wake_final, elems_coords_wake)
 
     def _getSectionsFromFile(self):
         """Obtain the nodes' location and row and cg of all elements.
@@ -780,10 +815,8 @@ class SolversInterface:
                     for i in range(len(cgSec)):
                         cgSec[i,:3] = (nodesSec[i,:3] + nodesSec[i+1,:3])/2
                         cgSec[i, 3] = i
-
-                    self.vBnd[ibody][iy][iReg].initStructures(nodesSec.shape[0], cgSec.shape[0])
-                    self.vBnd[ibody][iy][iReg].nodesCoord = nodesSec
-                    self.vBnd[ibody][iy][iReg].elemsCoord = cgSec
+                    
+                    self.vBnd[ibody][iy][iReg].setMesh(nodesSec[:,:3], cgSec[:,:3])
 
     def _checkInput(self, cfg:dict)->dict:
         """Check the input configuration.
@@ -840,9 +873,9 @@ class SolversInterface:
                 else:
                     iname = 'iWing'
                     vname = 'vAirfoil'
-                iBnd[ibody].append(Group(iname, typ, self.getnDim()))
+                iBnd[ibody].append(InviscidData(iname, typ, self.getnDim()))
                 for isec, s in enumerate(self.cfg['sections'][bodyName]):
-                        vBnd[ibody][isec].append(Group(vname, typ, self.getnDim()))
+                        vBnd[ibody][isec].append(ViscousData(vname, typ, self.getnDim()))
         return iBnd, vBnd
 
     def _getInterpolator(self)->object:
@@ -860,11 +893,13 @@ class SolversInterface:
                 for iSec in range(len(self.vBnd)):
                     for iReg, vreg in enumerate(self.vBnd[ibody][iSec]):
                         ireg = self.iBnd[ibody][iReg]
-                        vreg.initStructures(ireg.nPoints, ireg.nElems)
-                        vreg.nodesCoord = ireg.nodesCoord.copy()
-                        vreg.elemsCoord = ireg.elemsCoord.copy()
-                        vreg.connectListNodes = ireg.connectListNodes.copy()
-                        vreg.connectListElems = ireg.connectListElems.copy()
+                        nodes = ireg.getNodes('all').copy()
+                        elems = ireg.getElms('all').copy()
+                        vreg.setMesh(nodes, elems)
+                        vreg.setConnectLists(ireg._connectNodes.copy(), ireg._connectElems.copy())
+                        vreg.setRows(ireg.getNodeRows(), ireg.getElmRows())
+                        # vreg.connectListNodes = ireg._connectNodes.copy()
+                        # vreg.connectListElems = ireg._connectElems.copy()
 
         elif self.cfg['interpolator'] == 'Rbf':
             from blast.interfaces.interpolators.blRbfInterpolator import RbfInterpolator as interp
@@ -910,20 +945,24 @@ class SolversInterface:
             # Get all named regions
             allregs = []
             for vreg in self.vBnd[ibody][0]:
-                allregs.append(vreg.type)
+                allregs.append(vreg.getType())
             # Compute average and max AR of viscous cells
             avgAR = 0
             nNodes = 0
             if self.getnDim() == 3:
                 point_cpt = 0
                 for isec in range(1,len(self.vBnd[ibody])):
-                    if self.vBnd[ibody][isec][0].nodesCoord.shape[0] != self.vBnd[ibody][isec-1][0].nodesCoord.shape[0]:
-                        # Sections do not have the same number of elements
+                    wing_section = self.vBnd[ibody][isec][0]
+                    wing_section_prev = self.vBnd[ibody][isec-1][0]
+                    if wing_section.getnNodes() != wing_section_prev.getnNodes():
+                        # Sections do not have the same number of nodes
                         point_cpt = 1
                         continue
-                    for ipoint in range(1,self.vBnd[ibody][isec][0].nodesCoord.shape[0]):
-                        dx = abs((self.vBnd[ibody][isec][0].nodesCoord[ipoint,0] - self.vBnd[ibody][isec][0].nodesCoord[ipoint-1,0]))
-                        dy = abs((self.vBnd[ibody][isec][0].nodesCoord[ipoint,1] - self.vBnd[ibody][isec-1][0].nodesCoord[ipoint,1]))
+                    nodes_sec = wing_section.getNodes('all')
+                    nodes_prev = wing_section_prev.getNodes('all')
+                    for ipoint in range(1,wing_section.getnNodes()):
+                        dx = abs((nodes_sec[ipoint,0] - nodes_sec[ipoint-1,0]))
+                        dy = abs((nodes_sec[ipoint,1] - nodes_prev[ipoint,1]))
                         # AR = b^2 / S
                         avgAR += dy**2 / (dx * dy)
                         point_cpt += 1
@@ -931,10 +970,86 @@ class SolversInterface:
                 avgAR /= point_cpt
             elif self.getnDim() == 2:
                 for isec in range(len(self.vBnd[ibody])):
-                    for ipoint in range(1,self.vBnd[ibody][isec][0].nodesCoord.shape[0]):
+                    for ipoint in range(1,self.vBnd[ibody][isec][0].getnNodes()):
                         nNodes += 1
             print(f'  - {bodyName}: {len(self.vBnd[ibody])} section(s) with sides {allregs},\n\t{nNodes} surface nodes, avg AR {avgAR:.2f}')
         print('')
+    
+    def save_restart(self, cd, sfx=''):
+        import h5py
+
+        if 'restart_outputName' in self.cfg:
+            filename = self.cfg['restart_outputName']
+        else:
+            filename = 'blaster_restart'
+
+        if self.getVerbose() > 0:
+            print(f'Saving restart file {filename+sfx}.h5', end='...')
+        with h5py.File(f'{filename+sfx}.h5', 'w') as f:
+            f.attrs["cd"] = cd
+            f.attrs["nDim"] = self.getnDim()
+            for ibody in range(len(self.vBnd)):
+                body = f.create_group(f"body{ibody}")
+                for isec in range(len(self.vBnd[ibody])):
+                    section = body.create_group(f"section{isec}")
+                    for ireg in range(len(self.vBnd[ibody][isec])):
+                        pyreg = self.vBnd[ibody][isec][ireg]
+                        region = section.create_group(f"region{ireg}")
+                        region.attrs["name"] = pyreg.getName()
+                        region.attrs["type"] = pyreg.getType()
+                        if pyreg.getStagPoint() is None:
+                            region.attrs["stagnation_index"] = -1
+                        else:
+                            region.attrs["stagnation_index"] = pyreg.getStagPoint()
+                        region.create_dataset("nodes", data=pyreg.getNodes('all'))
+                        region.create_dataset("elems", data=pyreg.getElms('all'))
+                        region.create_dataset("blowing", data=pyreg.getBlowingVelocity('all'))
+                    for ireg in range(len(self.deltaStarExt[ibody][isec])):
+                        region = section.create_group(f"blInteraction{ireg}")
+                        region.create_dataset("deltaStarExt", data=self.deltaStarExt[ibody][isec][ireg])
+                        region.create_dataset("xxExt", data=self.xxExt[ibody][isec][ireg])
+                        region.create_dataset("vtExt", data=self.vtExt[ibody][isec][ireg])
+        if self.getVerbose() > 0:
+            print('done')
+    
+    def _load_restart(self, restart_file, vBnd):
+        # Override viscous initialization flag to avoid reinitializing the interaction laws
+        self.vinit = True
+        self._initializeInteractionLaw()
+
+        import h5py
+        with h5py.File(restart_file+'.h5', 'r') as f:
+            self.cdPrev = f.attrs["cd"]
+            if f.attrs["nDim"] != self.getnDim():
+                raise RuntimeError("Dimension mismatch in restart file")
+            for ibody in range(len(vBnd)):
+                body = f[f"body{ibody}"]
+                for isec in range(len(vBnd[ibody])):
+                    section = body[f"section{isec}"]
+                    # The regions are defined on wing and wake (no upper/lower separation)
+                    for ireg in range(len(vBnd[ibody][isec])):
+                        pyreg = vBnd[ibody][isec][ireg]
+                        region = section[f"region{ireg}"]
+                        if pyreg.getName() != region.attrs["name"]:
+                            raise RuntimeError(f"Name mismatch in restart file for body {ibody}, section {isec}, region {ireg}, expected {pyreg.getName()}, got {region['name'][()]}")
+                        if pyreg.getType() != region.attrs["type"]:
+                            raise RuntimeError(f"Type mismatch in restart file for body {ibody}, section {isec}, region {ireg}, expected {pyreg.getType()}, got {region['type'][()]}")
+                        if region.attrs["stagnation_index"] == -1:
+                            stag = None # Wake
+                        else:
+                            stag = region.attrs["stagnation_index"]
+                        pyreg.setStagPoint(stag)
+                        nodes = region["nodes"][:]
+                        elems = region["elems"][:]
+                        pyreg.setMesh(nodes, elems)
+                        pyreg.setBlowingVelocity(region["blowing"])
+                    # Coupling parameters are defined on upper/lower sides and wake
+
+            for iReg in range(len(self.deltaStarExt[ibody][isec])):
+                blregion = section[f"blInteraction{iReg}"]
+                self.deltaStarExt[ibody][isec][iReg] = blregion["deltaStarExt"][:]
+                self.vtExt[ibody][isec][iReg] = blregion["vtExt"][:]
+                self.xxExt[ibody][isec][iReg] = blregion["xxExt"][:]
 
     # Abstract methods for the interface
     def run()->int:
diff --git a/blast/tests/t_restart_2D.py b/blast/tests/t_restart_2D.py
new file mode 100644
index 0000000000000000000000000000000000000000..d4df7880f364535e150d3ff5792a465f10857bd3
--- /dev/null
+++ b/blast/tests/t_restart_2D.py
@@ -0,0 +1,168 @@
+#!/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 2025
+# Test the restart file
+
+# Imports.
+
+import blast.blUtils as vutils
+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, 'growthRatio': 1.3, 'msTe' : 0.01, 'msLe' : 0.01}, # 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
+    'Upstream' : 'upstream',
+    # Freestream
+    'M_inf' : 0.2, # freestream Mach number
+    'AoA' : 1., # 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
+        'CFL0' : 1,                 # Inital CFL number of the calculation
+        'couplIter': 50,            # 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.
+        'xtrF' : [None, 0.4],       # Forced transition locations
+        'interpolator' : 'Matching',# Interpolator for the coupling
+        
+        # Restart file
+        'restart_solution' : False, # Restart solution from file
+        'restart_inputName': None, # Name of the restart file to read
+        'restart_outputName': 'naca0012_restart', # Name of the restart file saved on disk
+        'iterSaveRestart' : 1, # Number of iterations between restart file save
+        'saveCouplingIters': True,
+    }
+
+def main():
+    import os.path
+
+    # Iteration to restart
+    iterRestart = 10
+    
+    # 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 = vutils.initBlast(icfg, vcfg)
+    tms['pre'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    aeroCoeffs_ref = coupler.run()
+    tms['solver'].stop()
+
+    del isol, vsol, coupler
+
+    vcfg["restart_solution"] = True
+    vcfg["restart_inputName"] = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))\
+                                + '/workspace/blast_tests_t_restart/'\
+                                + vcfg["restart_outputName"] + f"_{iterRestart}"
+    vcfg["iterSaveRestart"] = float('inf')
+
+
+    coupler, isol, vsol = vutils.initBlast(icfg, vcfg)
+    aeroCoeff_restart = coupler.run()
+
+    # Save in vtk format
+    vutils.save(vsol.bodies[0].sections)
+
+    # 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()))
+    print('')
+    total_iters_ref = len(aeroCoeffs_ref['Cl'])
+    total_iters_restart = len(aeroCoeff_restart['Cl'])
+    print(f"Number of iterations without restart: {total_iters_ref}")
+    print(f'Restarted at iteration {iterRestart}')
+    print(f"Number of iterations with restart: {total_iters_restart}")
+    tms['total'].stop()
+
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    print('SOLVERS statistics')
+    print(coupler.tms)
+
+    bl_solution = vutils.getSolution(vsol.bodies[0].sections, write=False)[0]
+
+    final_res_ref = np.log10((aeroCoeffs_ref['Cd'][-1] - aeroCoeffs_ref['Cd'][-2]) / aeroCoeffs_ref['Cd'][-1])
+    final_res_restart = np.log10((aeroCoeff_restart['Cd'][-1] - aeroCoeff_restart['Cd'][-2]) / aeroCoeff_restart['Cd'][-1])
+    
+    # Test solution
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    tests.add(CTest('Cl compare', aeroCoeff_restart['Cl'][-1], aeroCoeffs_ref['Cl'][-1], 1e-8))
+    tests.add(CTest('Cd compare', aeroCoeff_restart['Cd'][-1], aeroCoeffs_ref['Cd'][-1], 1e-8))
+    tests.add(CTest('nIters', total_iters_restart, total_iters_ref-iterRestart-1, 0, forceabs=True))
+    tests.add(CTest('Final residual Cd (log)', final_res_restart, final_res_ref, 1e-6))
+    tests.run()
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()