(feat) Modify section generation and interpolation

The upper and the lower side are now interpolated in two groups

blast/interfaces/DSolversInterface.py

@@ -31,12 +31,16 @@ class SolversInterface:
         elif self.cfg['interpolator'] == 'Rbf':
             from blast.interfaces.interpolators.DRbfInterpolator import RbfInterpolator as interp
             print('Initializing RBF interpolator.')
-            self.getVWing()
+            upNodes, lwNodes, upElms, lwElms = self.getVWing2()
         else:
             raise RuntimeError('Incorrect interpolator specified in config.')
 
         self.interpolator = interp(self.cfg, self.getnDim())
 
+        if self.cfg['interpolator'] == 'Rbf':
+            self.interpolator.setSides(upNodes, lwNodes, upElms, lwElms)
+            print('Setting sides in RBF interpolator.')
+
         # Initialize transfer quantities
         self.deltaStarExt = [[np.zeros(0) for iReg in range(3)]\
                              for iSec in range(self.cfg['nSections'])]
@@ -298,6 +302,171 @@ class SolversInterface:
                 self.vtExt[iSec][i] = np.zeros(reg.getnNodes())
         self.vinit = True
         return 0
+    
+    def getVWing2(self):
+        from scipy.interpolate import griddata
+        from scipy.interpolate import interp1d
+        import matplotlib.pyplot as plt
+
+        upper_elems = []
+        lower_elems = []
+        upper_nodes = np.zeros((0, 4))
+        lower_nodes = np.zeros((0, 4))
+        try:
+            lower_elems = np.loadtxt(self.blw[0].tag.elems[0].ptag.name + '_.dat', dtype=int)
+        except:
+            raise RuntimeError('Could not find file containing lower side of wing.')
+
+        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]:
+                        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]:
+                        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_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 upper_nodes.shape[0] == 0 or lower_nodes.shape[0] == 0:
+            raise RuntimeError('Could not identify upper or lower side.')
+
+        sections = self.cfg['sections']
+        interpolated_upper = []
+        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"""
+        
+        # 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)
+
+            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))
+
+            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
+        
+        # Wake
+        wake = np.zeros((0,3))
+        for n in self.blw[1].nodes:
+            wake = np.row_stack((wake, [n.pos[0], n.pos[1], n.pos[2]]))
+
+        sections = self.cfg['sections']
+        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
+
+        # 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)
+
+            wf = np.column_stack((x, np.full_like(x, sections[i]), z))
+            wake_final.append(np.flip(wf, axis=0))
+        
+        elems_coords_wake = []
+        for isec, sec in enumerate(wake_final):
+            elems_coords_wake.append(np.zeros((sec.shape[0]-1, 3)))
+            for i in range(sec.shape[0]-1):
+                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].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].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.

blast/interfaces/dart/DartInterface.py

@@ -363,7 +363,7 @@ class DartInterface(SolversInterface):
         """
 
         for iRegion in range(len(self.cfg['iMsh'])):
-            for e in self.cfg['iMsh'][iRegion].tag.elems:
+            for e in self.blw[iRegion].tag.elems:
                 # Compute cg position
                 cg_pt = np.zeros(3)
                 for nw in e.nodes:
@@ -376,16 +376,15 @@ 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))
+        # 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])

blast/interfaces/interpolators/DRbfInterpolator.py

@@ -7,7 +7,137 @@ class RbfInterpolator:
         self.cfg = _cfg
         self.nDim = _nDim
 
+    def setSides(self, _upNodes, _lwNodes, _upElems, _lwElems):
+        self.upNodes = _upNodes
+        self.lwNodes = _lwNodes
+        self.upElems = _upElems
+        self.lwElems = _lwElems
+
+
     def inviscidToViscous(self, iDict, vDict):
+        
+        print('Inviscid to viscous')
+        # upper side
+        x0_up = np.zeros((0, self.nDim+1))
+        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])
+        for i, row in enumerate(x0_up[:,3]):
+            M0_up[i] = iDict[0].M[iDict[0].nodesCoord[:,3] == row][0]
+        rho0_up = np.zeros(x0_up.shape[0])
+        for i, row in enumerate(x0_up[:,3]):
+            rho0_up[i] = iDict[0].Rho[iDict[0].nodesCoord[:,3] == row][0]
+        v0_up = np.zeros((x0_up.shape[0], self.nDim))
+        for i, row in enumerate(x0_up[:,3]):
+            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))
+        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])
+        for i, row in enumerate(x0_lw[:,3]):
+            M0_lw[i] = iDict[0].M[iDict[0].nodesCoord[:,3] == row][0]
+        rho0_lw = np.zeros(x0_lw.shape[0])
+        for i, row in enumerate(x0_lw[:,3]):
+            rho0_lw[i] = iDict[0].Rho[iDict[0].nodesCoord[:,3] == row][0]
+        v0_lw = np.zeros((x0_lw.shape[0], self.nDim))
+        for i, row in enumerate(x0_lw[:,3]):
+            for idim in range(self.nDim):
+                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)
+            v_up = np.zeros((xInterp_up.shape[0], self.nDim))
+            for idim in range(self.nDim):
+                v_up[:,idim] = self.__rbfToSection(x0_up[:,:self.nDim], v0_up[:,idim], xInterp_up)
+
+            xInterp_lw = vDict[iSec][0].nodesCoord[np.argmin(vDict[iSec][0].nodesCoord[:,0]):, :self.nDim]
+            M_lw = self.__rbfToSection(x0_lw[:,:self.nDim], M0_lw, xInterp_lw)
+            rho_lw = self.__rbfToSection(x0_lw[:,:self.nDim], rho0_lw, xInterp_lw)
+            v_lw = np.zeros((xInterp_lw.shape[0], self.nDim))
+            for idim in range(self.nDim):
+                v_lw[:,idim] = self.__rbfToSection(x0_lw[:,:self.nDim], v0_lw[:,idim], xInterp_lw)
+
+            xInterp_wk = vDict[iSec][1].nodesCoord[:, :self.nDim]
+            M_wk = self.__rbfToSection(iDict[1].nodesCoord[:,:self.nDim], iDict[1].M, xInterp_wk)
+            rho_wk = self.__rbfToSection(iDict[1].nodesCoord[:,:self.nDim], iDict[1].Rho, xInterp_wk)
+            v_wk = np.zeros((xInterp_wk.shape[0], self.nDim))
+            for idim in range(self.nDim):
+                v_wk[:,idim] = self.__rbfToSection(iDict[1].nodesCoord[:,:self.nDim], iDict[1].V[:,idim], xInterp_wk)
+            
+            M_int = np.concatenate((M_up, M_lw))
+            rho_int = np.concatenate((rho_up, rho_lw))
+            v_int = np.row_stack((v_up, v_lw))
+
+            vDict[iSec][0].updateVars(v_int, M_int, rho_int)
+
+            vDict[iSec][1].updateVars(v_wk, M_wk, rho_wk)
+            
+
+            # 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))
+            # for idim in range(self.nDim):
+            #     V_tot[:,idim] = self.__rbfToSection(iDict[0].nodesCoord[:,:self.nDim], iDict[0].V[:, idim], vDict[iSec][0].nodesCoord[:,:self.nDim])
+            # Rho_tot = self.__rbfToSection(iDict[0].nodesCoord[:,:self.nDim], iDict[0].Rho, vDict[iSec][0].nodesCoord[:,:self.nDim])
+
+            # from mpl_toolkits.mplot3d import Axes3D
+            # import matplotlib.pyplot as plt
+            # fig = plt.figure()
+            # ax = fig.add_subplot(111, projection='3d')
+            # ax.plot_trisurf(x0_up[:,0], x0_up[:,1], M0_up)
+            # ax.plot(xInterp_up[:,0], xInterp_up[:,1], M_up, lw = 2)
+            # plt.draw()
+            # fig = plt.figure()
+            # ax = fig.add_subplot(111, projection='3d')
+            # ax.plot_trisurf(x0_lw[:,0], x0_lw[:,1], M0_lw)
+            # ax.plot(xInterp_lw[:,0], xInterp_lw[:,1], M_lw, lw = 2)
+            # plt.draw()
+
+            # plt.figure()
+            # plt.plot(xInterp_up[:,0], M_up, label='upper')
+            # plt.plot(xInterp_lw[:,0], M_lw, label='lower')
+            # plt.plot(vDict[iSec][0].nodesCoord[:,0], M_tot, '--', label='total')
+            # plt.legend(frameon=False)
+            # plt.title('Mach number')
+            # plt.draw()
+
+            # plt.figure()
+            # plt.plot(xInterp_up[:,0], rho_up, label='upper')
+            # plt.plot(xInterp_lw[:,0], rho_lw, label='lower')
+            # plt.plot(vDict[iSec][0].nodesCoord[:,0], Rho_tot, '--', label='total')
+            # plt.legend(frameon=False)
+            # plt.title('Density')
+            # plt.draw()
+
+            # for idim in range(self.nDim):
+            #     plt.figure()
+            #     plt.plot(xInterp_up[:,0], v_up[:,idim], label='upper')
+            #     plt.plot(xInterp_lw[:,0], v_lw[:,idim], label='lower')
+            #     plt.plot(vDict[iSec][0].nodesCoord[:,0], V_tot[:,idim], '--', label='total')
+            #     plt.legend(frameon=False)
+            #     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)):
@@ -22,6 +152,58 @@ class RbfInterpolator:
                 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))
+        xv_wk = np.zeros((0, self.nDim))
+        blwv_up = np.zeros(0)
+        blwv_lw = np.zeros(0)
+        blwv_wk = np.zeros(0)
+        for iSec, sec in enumerate(vDict):
+            xv_up = np.row_stack((xv_up, sec[0].elemsCoord[:np.argmin(sec[0].nodesCoord[:,0])-1, :self.nDim]))
+            xv_lw = np.row_stack((xv_lw, sec[0].elemsCoord[np.argmin(sec[0].nodesCoord[:,0]):, :self.nDim]))
+            xv_wk = np.row_stack((xv_wk, sec[1].elemsCoord[:,:self.nDim]))
+            blwv_up = np.concatenate((blwv_up, sec[0].blowingVel[:np.argmin(sec[0].nodesCoord[:,0])-1]))
+            blwv_lw = np.concatenate((blwv_lw, sec[0].blowingVel[np.argmin(sec[0].nodesCoord[:,0]):]))
+            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))
+        for no in self.upElems:
+            xi_up = np.row_stack((xi_up, np.unique(iDict[0].elemsCoord[iDict[0].elemsCoord[:,3]==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_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))
+        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)])

blast/tests/dart/onera_3D.py

+#!/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()