diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index b7542bc56ab9b64d9600921f416f5607e88d3053..29e4081378ddbe88b673f6f82a8d1e312da5df70 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/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.
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index 7913ab549b1b745aa332d945590932b7f0633bca..ffa71953eba16204c3d6d41525f4603c1c7dd4b9 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/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])
diff --git a/blast/interfaces/interpolators/DRbfInterpolator.py b/blast/interfaces/interpolators/DRbfInterpolator.py
index 1eda4a3e12b525a09907cbdb0fede1241d7d51ba..4f5ae028116fb80d3fb6f0591660de1c0fe20816 100644
--- a/blast/interfaces/interpolators/DRbfInterpolator.py
+++ b/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)])
diff --git a/blast/models/dart/oneraM6.geo b/blast/models/dart/oneraM6.geo
new file mode 100644
index 0000000000000000000000000000000000000000..616835e51154fb213fbb6f5289e64be8f80a85c6
--- /dev/null
+++ b/blast/models/dart/oneraM6.geo
@@ -0,0 +1,591 @@
+/* Onera M6 wing */
+// Initially generated by unsRgridWingGen.m
+
+// Parameters
+// domain and mesh
+DefineConstant[ lgt = { 3.0, Name "Channel half length" } ];
+DefineConstant[ wdt = { 2.0, Name "Channel width" } ];
+DefineConstant[ hgt = { 3.0, Name "Channel half height" } ];
+DefineConstant[ msLeRt = { 0.01, Name "Root leading edge mesh size" } ];
+DefineConstant[ msTeRt = { 0.01, Name "Root trailing edge mesh size" } ];
+DefineConstant[ msLeTp = { 0.01, Name "Tip leading edge mesh size" } ];
+DefineConstant[ msTeTp = { 0.01, Name "Tip trailing edge mesh size" } ];
+DefineConstant[ msF = { 1.0, Name "Farfield mesh size" } ];
+
+
+//// GEOMETRY
+
+
+/// Points
+
+// Airfoil 1: oneraM6, 143 points
+Point(1) = {0.805900,0.000000,0.000000,msTeRt};
+Point(2) = {0.804129,0.000000,0.000232};
+Point(3) = {0.802038,0.000000,0.000505};
+Point(4) = {0.799569,0.000000,0.000828};
+Point(5) = {0.796652,0.000000,0.001210};
+Point(6) = {0.793208,0.000000,0.001661};
+Point(7) = {0.789139,0.000000,0.002193};
+Point(8) = {0.784331,0.000000,0.002822};
+Point(9) = {0.778649,0.000000,0.003565};
+Point(10) = {0.771932,0.000000,0.004444};
+Point(11) = {0.763989,0.000000,0.005483};
+Point(12) = {0.754592,0.000000,0.006710};
+Point(13) = {0.743470,0.000000,0.008151};
+Point(14) = {0.730299,0.000000,0.009828,msTeRt};
+Point(15) = {0.714691,0.000000,0.011690};
+Point(16) = {0.696182,0.000000,0.013774};
+Point(17) = {0.677546,0.000000,0.015783};
+Point(18) = {0.658809,0.000000,0.017717};
+Point(19) = {0.639970,0.000000,0.019586};
+Point(20) = {0.621026,0.000000,0.021397};
+Point(21) = {0.601980,0.000000,0.023159};
+Point(22) = {0.582826,0.000000,0.024875};
+Point(23) = {0.563568,0.000000,0.026547};
+Point(24) = {0.544202,0.000000,0.028170};
+Point(25) = {0.524729,0.000000,0.029737};
+Point(26) = {0.505147,0.000000,0.031238};
+Point(27) = {0.485455,0.000000,0.032658};
+Point(28) = {0.465652,0.000000,0.033984};
+Point(29) = {0.445738,0.000000,0.035197};
+Point(30) = {0.425711,0.000000,0.036282};
+Point(31) = {0.405570,0.000000,0.037225};
+Point(32) = {0.385317,0.000000,0.038011};
+Point(33) = {0.364947,0.000000,0.038631};
+Point(34) = {0.344460,0.000000,0.039073};
+Point(35) = {0.323856,0.000000,0.039344};
+Point(36) = {0.303135,0.000000,0.039432};
+Point(37) = {0.282293,0.000000,0.039343};
+Point(38) = {0.261331,0.000000,0.039078};
+Point(39) = {0.240248,0.000000,0.038642,1.5*msLeRt};
+Point(40) = {0.219042,0.000000,0.038037};
+Point(41) = {0.197712,0.000000,0.037261};
+Point(42) = {0.176258,0.000000,0.036306};
+Point(43) = {0.154679,0.000000,0.035154};
+Point(44) = {0.132972,0.000000,0.033774};
+Point(45) = {0.111131,0.000000,0.032117};
+Point(46) = {0.092662,0.000000,0.030457};
+Point(47) = {0.077073,0.000000,0.028830};
+Point(48) = {0.063937,0.000000,0.027269};
+Point(49) = {0.052891,0.000000,0.025800};
+Point(50) = {0.043619,0.000000,0.024441};
+Point(51) = {0.035851,0.000000,0.023203};
+Point(52) = {0.029356,0.000000,0.022027};
+Point(53) = {0.023936,0.000000,0.020812};
+Point(54) = {0.019428,0.000000,0.019503};
+Point(55) = {0.015684,0.000000,0.018096};
+Point(56) = {0.012586,0.000000,0.016619};
+Point(57) = {0.010032,0.000000,0.015114};
+Point(58) = {0.007931,0.000000,0.013618};
+Point(59) = {0.006214,0.000000,0.012165};
+Point(60) = {0.004815,0.000000,0.010776};
+Point(61) = {0.003683,0.000000,0.009463};
+Point(62) = {0.002775,0.000000,0.008233};
+Point(63) = {0.002050,0.000000,0.007089};
+Point(64) = {0.001480,0.000000,0.006027};
+Point(65) = {0.001037,0.000000,0.005045};
+Point(66) = {0.000698,0.000000,0.004138};
+Point(67) = {0.000444,0.000000,0.003301};
+Point(68) = {0.000260,0.000000,0.002529};
+Point(69) = {0.000135,0.000000,0.001818};
+Point(70) = {0.000056,0.000000,0.001162};
+Point(71) = {0.000013,0.000000,0.000557};
+Point(72) = {0.000000,0.000000,0.000000,msLeRt};
+Point(73) = {0.000013,0.000000,-0.000557};
+Point(74) = {0.000056,0.000000,-0.001162};
+Point(75) = {0.000135,0.000000,-0.001818};
+Point(76) = {0.000260,0.000000,-0.002529};
+Point(77) = {0.000444,0.000000,-0.003301};
+Point(78) = {0.000698,0.000000,-0.004138};
+Point(79) = {0.001037,0.000000,-0.005045};
+Point(80) = {0.001480,0.000000,-0.006027};
+Point(81) = {0.002050,0.000000,-0.007089};
+Point(82) = {0.002775,0.000000,-0.008233};
+Point(83) = {0.003683,0.000000,-0.009463};
+Point(84) = {0.004815,0.000000,-0.010776};
+Point(85) = {0.006214,0.000000,-0.012165};
+Point(86) = {0.007931,0.000000,-0.013618};
+Point(87) = {0.010032,0.000000,-0.015114};
+Point(88) = {0.012586,0.000000,-0.016619};
+Point(89) = {0.015684,0.000000,-0.018096};
+Point(90) = {0.019428,0.000000,-0.019503};
+Point(91) = {0.023936,0.000000,-0.020812};
+Point(92) = {0.029356,0.000000,-0.022027};
+Point(93) = {0.035851,0.000000,-0.023203};
+Point(94) = {0.043619,0.000000,-0.024441};
+Point(95) = {0.052891,0.000000,-0.025800};
+Point(96) = {0.063937,0.000000,-0.027269};
+Point(97) = {0.077073,0.000000,-0.028830};
+Point(98) = {0.092662,0.000000,-0.030457};
+Point(99) = {0.111131,0.000000,-0.032117};
+Point(100) = {0.132972,0.000000,-0.033774};
+Point(101) = {0.154679,0.000000,-0.035154};
+Point(102) = {0.176258,0.000000,-0.036306};
+Point(103) = {0.197712,0.000000,-0.037261};
+Point(104) = {0.219042,0.000000,-0.038037};
+Point(105) = {0.240248,0.000000,-0.038642,1.5*msLeRt};
+Point(106) = {0.261331,0.000000,-0.039078};
+Point(107) = {0.282293,0.000000,-0.039343};
+Point(108) = {0.303135,0.000000,-0.039432};
+Point(109) = {0.323856,0.000000,-0.039344};
+Point(110) = {0.344460,0.000000,-0.039073};
+Point(111) = {0.364947,0.000000,-0.038631};
+Point(112) = {0.385317,0.000000,-0.038011};
+Point(113) = {0.405570,0.000000,-0.037225};
+Point(114) = {0.425711,0.000000,-0.036282};
+Point(115) = {0.445738,0.000000,-0.035197};
+Point(116) = {0.465652,0.000000,-0.033984};
+Point(117) = {0.485455,0.000000,-0.032658};
+Point(118) = {0.505147,0.000000,-0.031238};
+Point(119) = {0.524729,0.000000,-0.029737};
+Point(120) = {0.544202,0.000000,-0.028170};
+Point(121) = {0.563568,0.000000,-0.026547};
+Point(122) = {0.582826,0.000000,-0.024875};
+Point(123) = {0.601980,0.000000,-0.023159};
+Point(124) = {0.621026,0.000000,-0.021397};
+Point(125) = {0.639970,0.000000,-0.019586};
+Point(126) = {0.658809,0.000000,-0.017717};
+Point(127) = {0.677546,0.000000,-0.015783};
+Point(128) = {0.696182,0.000000,-0.013774};
+Point(129) = {0.714691,0.000000,-0.011690};
+Point(130) = {0.730299,0.000000,-0.009828,msTeRt};
+Point(131) = {0.743470,0.000000,-0.008151};
+Point(132) = {0.754592,0.000000,-0.006710};
+Point(133) = {0.763989,0.000000,-0.005483};
+Point(134) = {0.771932,0.000000,-0.004444};
+Point(135) = {0.778649,0.000000,-0.003565};
+Point(136) = {0.784331,0.000000,-0.002822};
+Point(137) = {0.789139,0.000000,-0.002193};
+Point(138) = {0.793208,0.000000,-0.001661};
+Point(139) = {0.796652,0.000000,-0.001210};
+Point(140) = {0.799569,0.000000,-0.000828};
+Point(141) = {0.802038,0.000000,-0.000505};
+Point(142) = {0.804129,0.000000,-0.000232,msTeRt};
+
+// Airfoil 2: oneraM6, 143 points
+Point(144) = {1.143427,1.196000,0.000000,msTeTp};
+Point(145) = {1.142432,1.196000,0.000130};
+Point(146) = {1.141256,1.196000,0.000284};
+Point(147) = {1.139869,1.196000,0.000466};
+Point(148) = {1.138230,1.196000,0.000680};
+Point(149) = {1.136294,1.196000,0.000933};
+Point(150) = {1.134007,1.196000,0.001232};
+Point(151) = {1.131305,1.196000,0.001586};
+Point(152) = {1.128112,1.196000,0.002004};
+Point(153) = {1.124337,1.196000,0.002498};
+Point(154) = {1.119873,1.196000,0.003082};
+Point(155) = {1.114592,1.196000,0.003771};
+Point(156) = {1.108341,1.196000,0.004581};
+Point(157) = {1.100939,1.196000,0.005523,1.5*msTeTp};
+Point(158) = {1.092167,1.196000,0.006570};
+Point(159) = {1.081765,1.196000,0.007741};
+Point(160) = {1.071292,1.196000,0.008870};
+Point(161) = {1.060762,1.196000,0.009957};
+Point(162) = {1.050174,1.196000,0.011007};
+Point(163) = {1.039528,1.196000,0.012025};
+Point(164) = {1.028824,1.196000,0.013015};
+Point(165) = {1.018059,1.196000,0.013980};
+Point(166) = {1.007236,1.196000,0.014919};
+Point(167) = {0.996353,1.196000,0.015831};
+Point(168) = {0.985409,1.196000,0.016712};
+Point(169) = {0.974403,1.196000,0.017556};
+Point(170) = {0.963336,1.196000,0.018354};
+Point(171) = {0.952208,1.196000,0.019099};
+Point(172) = {0.941016,1.196000,0.019781};
+Point(173) = {0.929760,1.196000,0.020391};
+Point(174) = {0.918441,1.196000,0.020920};
+Point(175) = {0.907059,1.196000,0.021362};
+Point(176) = {0.895611,1.196000,0.021711};
+Point(177) = {0.884097,1.196000,0.021959};
+Point(178) = {0.872518,1.196000,0.022111};
+Point(179) = {0.860873,1.196000,0.022161};
+Point(180) = {0.849160,1.196000,0.022111};
+Point(181) = {0.837379,1.196000,0.021962};
+Point(182) = {0.825530,1.196000,0.021717,1.5*msLeTp};
+Point(183) = {0.813612,1.196000,0.021377};
+Point(184) = {0.801625,1.196000,0.020941};
+Point(185) = {0.789568,1.196000,0.020404};
+Point(186) = {0.777440,1.196000,0.019757};
+Point(187) = {0.765241,1.196000,0.018981};
+Point(188) = {0.752966,1.196000,0.018050};
+Point(189) = {0.742587,1.196000,0.017117};
+Point(190) = {0.733826,1.196000,0.016203};
+Point(191) = {0.726444,1.196000,0.015325};
+Point(192) = {0.720236,1.196000,0.014500};
+Point(193) = {0.715025,1.196000,0.013736};
+Point(194) = {0.710659,1.196000,0.013040};
+Point(195) = {0.707009,1.196000,0.012379};
+Point(196) = {0.703963,1.196000,0.011696};
+Point(197) = {0.701429,1.196000,0.010961};
+Point(198) = {0.699325,1.196000,0.010170};
+Point(199) = {0.697584,1.196000,0.009340};
+Point(200) = {0.696149,1.196000,0.008494};
+Point(201) = {0.694968,1.196000,0.007654};
+Point(202) = {0.694003,1.196000,0.006837};
+Point(203) = {0.693217,1.196000,0.006056};
+Point(204) = {0.692581,1.196000,0.005318};
+Point(205) = {0.692070,1.196000,0.004627};
+Point(206) = {0.691663,1.196000,0.003984};
+Point(207) = {0.691343,1.196000,0.003387};
+Point(208) = {0.691094,1.196000,0.002835};
+Point(209) = {0.690903,1.196000,0.002325};
+Point(210) = {0.690760,1.196000,0.001855};
+Point(211) = {0.690657,1.196000,0.001421};
+Point(212) = {0.690587,1.196000,0.001022};
+Point(213) = {0.690542,1.196000,0.000653};
+Point(214) = {0.690518,1.196000,0.000313};
+Point(215) = {0.690511,1.196000,0.000000,msLeTp};
+Point(216) = {0.690518,1.196000,-0.000313};
+Point(217) = {0.690542,1.196000,-0.000653};
+Point(218) = {0.690587,1.196000,-0.001022};
+Point(219) = {0.690657,1.196000,-0.001421};
+Point(220) = {0.690760,1.196000,-0.001855};
+Point(221) = {0.690903,1.196000,-0.002325};
+Point(222) = {0.691094,1.196000,-0.002835};
+Point(223) = {0.691343,1.196000,-0.003387};
+Point(224) = {0.691663,1.196000,-0.003984};
+Point(225) = {0.692070,1.196000,-0.004627};
+Point(226) = {0.692581,1.196000,-0.005318};
+Point(227) = {0.693217,1.196000,-0.006056};
+Point(228) = {0.694003,1.196000,-0.006837};
+Point(229) = {0.694968,1.196000,-0.007654};
+Point(230) = {0.696149,1.196000,-0.008494};
+Point(231) = {0.697584,1.196000,-0.009340};
+Point(232) = {0.699325,1.196000,-0.010170};
+Point(233) = {0.701429,1.196000,-0.010961};
+Point(234) = {0.703963,1.196000,-0.011696};
+Point(235) = {0.707009,1.196000,-0.012379};
+Point(236) = {0.710659,1.196000,-0.013040};
+Point(237) = {0.715025,1.196000,-0.013736};
+Point(238) = {0.720236,1.196000,-0.014500};
+Point(239) = {0.726444,1.196000,-0.015325};
+Point(240) = {0.733826,1.196000,-0.016203};
+Point(241) = {0.742587,1.196000,-0.017117};
+Point(242) = {0.752966,1.196000,-0.018050};
+Point(243) = {0.765241,1.196000,-0.018981};
+Point(244) = {0.777440,1.196000,-0.019757};
+Point(245) = {0.789568,1.196000,-0.020404};
+Point(246) = {0.801625,1.196000,-0.020941};
+Point(247) = {0.813612,1.196000,-0.021377};
+Point(248) = {0.825530,1.196000,-0.021717,1.5*msLeTp};
+Point(249) = {0.837379,1.196000,-0.021962};
+Point(250) = {0.849160,1.196000,-0.022111};
+Point(251) = {0.860873,1.196000,-0.022161};
+Point(252) = {0.872518,1.196000,-0.022111};
+Point(253) = {0.884097,1.196000,-0.021959};
+Point(254) = {0.895611,1.196000,-0.021711};
+Point(255) = {0.907059,1.196000,-0.021362};
+Point(256) = {0.918441,1.196000,-0.020920};
+Point(257) = {0.929760,1.196000,-0.020391};
+Point(258) = {0.941016,1.196000,-0.019781};
+Point(259) = {0.952208,1.196000,-0.019099};
+Point(260) = {0.963336,1.196000,-0.018354};
+Point(261) = {0.974403,1.196000,-0.017556};
+Point(262) = {0.985409,1.196000,-0.016712};
+Point(263) = {0.996353,1.196000,-0.015831};
+Point(264) = {1.007236,1.196000,-0.014919};
+Point(265) = {1.018059,1.196000,-0.013980};
+Point(266) = {1.028824,1.196000,-0.013015};
+Point(267) = {1.039528,1.196000,-0.012025};
+Point(268) = {1.050174,1.196000,-0.011007};
+Point(269) = {1.060762,1.196000,-0.009957};
+Point(270) = {1.071292,1.196000,-0.008870};
+Point(271) = {1.081765,1.196000,-0.007741};
+Point(272) = {1.092167,1.196000,-0.006570};
+Point(273) = {1.100939,1.196000,-0.005523,1.5*msTeTp};
+Point(274) = {1.108341,1.196000,-0.004581};
+Point(275) = {1.114592,1.196000,-0.003771};
+Point(276) = {1.119873,1.196000,-0.003082};
+Point(277) = {1.124337,1.196000,-0.002498};
+Point(278) = {1.128112,1.196000,-0.002004};
+Point(279) = {1.131305,1.196000,-0.001586};
+Point(280) = {1.134007,1.196000,-0.001232};
+Point(281) = {1.136294,1.196000,-0.000933};
+Point(282) = {1.138230,1.196000,-0.000680};
+Point(283) = {1.139869,1.196000,-0.000466};
+Point(284) = {1.141256,1.196000,-0.000284};
+Point(285) = {1.142432,1.196000,-0.000130,msTeTp};
+
+// Box:
+Point(5001) = {-lgt,0.000000,-hgt,msF};
+Point(5002) = {1+lgt,0.000000,-hgt,msF};
+Point(5003) = {-lgt,0.000000,hgt,msF};
+Point(5004) = {1+lgt,0.000000,hgt,msF};
+Point(5005) = {-lgt,wdt,-hgt,msF};
+Point(5006) = {1+lgt,wdt,-hgt,msF};
+Point(5007) = {-lgt,wdt,hgt,msF};
+Point(5008) = {1+lgt,wdt,hgt,msF};
+
+// Tip:
+Point(5101) = {1.142432,1.196000,0.000000};
+Point(5102) = {1.141256,1.196040,0.000000};
+Point(5103) = {1.139869,1.196088,0.000000};
+Point(5104) = {1.138230,1.196144,0.000000};
+Point(5105) = {1.136294,1.196210,0.000000};
+Point(5106) = {1.134007,1.196289,0.000000};
+Point(5107) = {1.131305,1.196381,0.000000};
+Point(5108) = {1.128112,1.196491,0.000000};
+Point(5109) = {1.124337,1.196620,0.000000};
+Point(5110) = {1.119873,1.196773,0.000000};
+Point(5111) = {1.114592,1.196954,0.000000};
+Point(5112) = {1.108341,1.197168,0.000000};
+Point(5113) = {1.100939,1.197422,0.000000,1.5*msTeTp};
+Point(5114) = {1.092167,1.197722,0.000000};
+Point(5115) = {1.081765,1.198079,0.000000};
+Point(5116) = {1.071292,1.198438,0.000000};
+Point(5117) = {1.060762,1.198798,0.000000};
+Point(5118) = {1.050174,1.199161,0.000000};
+Point(5119) = {1.039528,1.199526,0.000000};
+Point(5120) = {1.028824,1.199893,0.000000};
+Point(5121) = {1.018059,1.200262,0.000000};
+Point(5122) = {1.007236,1.200632,0.000000};
+Point(5123) = {0.996353,1.201005,0.000000};
+Point(5124) = {0.985409,1.201380,0.000000};
+Point(5125) = {0.974403,1.201757,0.000000};
+Point(5126) = {0.963336,1.202137,0.000000};
+Point(5127) = {0.952208,1.202518,0.000000};
+Point(5128) = {0.941016,1.202901,0.000000};
+Point(5129) = {0.929760,1.203287,0.000000};
+Point(5130) = {0.918441,1.203675,0.000000};
+Point(5131) = {0.907059,1.204065,0.000000};
+Point(5132) = {0.895611,1.204457,0.000000};
+Point(5133) = {0.884097,1.204852,0.000000};
+Point(5134) = {0.872518,1.205248,0.000000};
+Point(5135) = {0.860873,1.205648,0.000000};
+Point(5136) = {0.849160,1.206049,0.000000};
+Point(5137) = {0.837379,1.206453,0.000000};
+Point(5138) = {0.825530,1.206859,0.000000,1.5*msLeTp};
+Point(5139) = {0.699294,1.211213,0.000000};
+
+// Dummy tip center:
+Point(5349) = {1.100939,1.196000,0.000000};
+Point(5350) = {0.825530,1.196000,0.000000};
+
+
+// Midplane:
+Point(5351) = {1+lgt,0.000000,0.000000,msF};
+Point(5352) = {1+lgt,1.196000,0.000000,msF};
+Point(5353) = {1+lgt,wdt,0.000000,msF};
+Point(5354) = {1.143427,wdt,0.000000,msF};
+Point(5355) = {1.100939,wdt,0.000000,msF};
+Point(5356) = {0.825530,wdt,0.000000,msF};
+Point(5357) = {0.690511,wdt,0.000000,msF};
+Point(5358) = {-lgt,wdt,0.000000,msF};
+Point(5359) = {-lgt,1.196000,0.000000,msF};
+Point(5360) = {-lgt,0.000000,0.000000,msF};
+
+/// Lines
+
+// Airfoil 1:
+Spline(1) = {1,2,3,4,5,6,7,8,9,10,11,12,13,14};
+Spline(2) = {14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39};
+Spline(3) = {39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72};
+Spline(4) = {72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105};
+Spline(5) = {105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130};
+Spline(6) = {130,131,132,133,134,135,136,137,138,139,140,141,142,1};
+
+// Airfoil 2:
+Spline(7) = {144,145,146,147,148,149,150,151,152,153,154,155,156,157};
+Spline(8) = {157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182};
+Spline(9) = {182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215};
+Spline(10) = {215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248};
+Spline(11) = {248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273};
+Spline(12) = {273,274,275,276,277,278,279,280,281,282,283,284,285,144};
+
+
+// Box:
+Line(41) = {5001,5002};
+Line(42) = {5003,5004};
+Line(43) = {5005,5006};
+Line(44) = {5007,5008};
+Line(45) = {5001,5005};
+Line(46) = {5002,5006};
+Line(47) = {5003,5007};
+Line(48) = {5004,5008};
+Line(49) = {5001,5360};
+Line(50) = {5002,5351};
+Line(51) = {5003,5360};
+Line(52) = {5004,5351};
+Line(53) = {5005,5358};
+Line(54) = {5006,5353};
+Line(55) = {5007,5358};
+Line(56) = {5008,5353};
+
+// Airfoil 1 to airfoil 2:
+Line(61) = {1,144};
+Line(62) = {14,157};
+Line(63) = {39,182};
+Line(64) = {72,215};
+Line(65) = {105,248};
+Line(66) = {130,273};
+
+
+// Tip:
+Spline(121) = {144,5101,5102,5103,5104,5105,5106,5107,5108,5109,5110,5111,5112,5113};
+Spline(122) = {5113,5114,5115,5116,5117,5118,5119,5120,5121,5122,5123,5124,5125,5126,5127,5128,5129,5130,5131,5132,5133,5134,5135,5136,5137,5138};
+If(GMSH_MAJOR_VERSION >= 4)
+ Bezier(123) = {5138,5139,215};
+Else
+ BSpline(123) = {5138,5139,215};
+EndIf
+Ellipse(124) = {157,5349,157,5113};
+Ellipse(125) = {182,5350,182,5138};
+Ellipse(126) = {248,5350,248,5138};
+Ellipse(127) = {273,5349,273,5113};
+
+// Midplane:
+Line(131) = {5351,5352};
+Line(132) = {5352,5353};
+Line(133) = {5353,5354};
+Line(134) = {5354,5355};
+Line(135) = {5355,5356};
+Line(136) = {5356,5357};
+Line(137) = {5357,5358};
+Line(138) = {5358,5359};
+Line(139) = {5359,5360};
+
+// Wing to midplane:
+Line(161) = {1,5351};
+Line(162) = {144,5352};
+Line(163) = {144,5354};
+Line(164) = {5113,5355};
+Line(165) = {5138,5356};
+Line(166) = {215,5357};
+Line(167) = {215,5359};
+Line(168) = {72,5360};
+
+/// Line loops & Surfaces
+
+// Box:
+Line Loop(1) = {1,2,3,168,-51,42,52,-161};
+Line Loop(2) = {48,56,-132,-131,-52};
+Line Loop(3) = {44,56,133,134,135,136,137,-55};
+Line Loop(4) = {47,55,138,139,-51};
+Line Loop(5) = {42,48,-44,-47};
+Line Loop(6) = {-6,-5,-4,168,-49,41,50,-161};
+Line Loop(7) = {46,54,-132,-131,-50};
+Line Loop(8) = {43,54,133,134,135,136,137,-53};
+Line Loop(9) = {45,53,138,139,-49};
+Line Loop(10) = {41,46,-43,-45};
+Plane Surface(1) = {1};
+Plane Surface(2) = {2};
+Plane Surface(3) = {-3};
+Plane Surface(4) = {-4};
+Plane Surface(5) = {-5};
+Plane Surface(6) = {-6};
+Plane Surface(7) = {-7};
+Plane Surface(8) = {8};
+Plane Surface(9) = {9};
+Plane Surface(10) = {10};
+
+// Wing 1:
+Line Loop(11) = {1,62,-7,-61};
+Line Loop(12) = {2,63,-8,-62};
+Line Loop(13) = {3,64,-9,-63};
+Line Loop(14) = {4,65,-10,-64};
+Line Loop(15) = {5,66,-11,-65};
+Line Loop(16) = {6,61,-12,-66};
+Surface(11) = {-11};
+Surface(12) = {-12};
+Surface(13) = {-13};
+Surface(14) = {-14};
+Surface(15) = {-15};
+Surface(16) = {-16};
+
+// Wingtip:
+Line Loop(71) = {7,124,-121};
+Line Loop(72) = {8,125,-122,-124};
+Line Loop(73) = {9,-123,-125};
+Line Loop(74) = {10,126,123};
+Line Loop(75) = {11,127,122,-126};
+Line Loop(76) = {12,121,-127};
+Surface(71) = {-71};
+Surface(72) = {-72};
+Surface(73) = {-73};
+Surface(74) = {-74};
+Surface(75) = {-75};
+Surface(76) = {-76};
+
+// Midplane:
+Line Loop(81) = {161,131,-162,-61};
+Line Loop(82) = {162,132,133,-163};
+Line Loop(83) = {163,134,-164,-121};
+Line Loop(84) = {164,135,-165,-122};
+Line Loop(85) = {165,136,-166,-123};
+Line Loop(86) = {167,-138,-137,-166};
+Line Loop(87) = {167,139,-168,64};
+Surface(81) = {81};
+Surface(82) = {82};
+Surface(83) = {83};
+Surface(84) = {84};
+Surface(85) = {85};
+Surface(86) = {-86};
+Surface(87) = {87};
+
+/// Surface loops & Volumes
+
+// Upper:
+Surface Loop(1) = {11,12,13,71,72,73,1,2,3,4,5,81,82,83,84,85,86,87};
+Volume(1) = {1};
+// Lower:
+Surface Loop(2) = {14,15,16,74,75,76,6,7,8,9,10,81,82,83,84,85,86,87};
+Volume(2) = {2};
+
+
+
+//// MESHING ALGORITHM
+
+
+/// 2D:
+
+///Wing, farfield and symmetry plane:
+MeshAlgorithm Surface{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,81,82,83,84,85,86,87} = 5;
+
+///Tip:
+MeshAlgorithm Surface{73,74} = 1;
+MeshAlgorithm Surface{71,72,75,76} = 5;
+
+/// 3D:
+
+Mesh.Algorithm3D = 2;
+Mesh.Optimize = 1;
+Mesh.Smoothing = 10;
+Mesh.SmoothNormals = 1;
+
+
+
+//// PHYSICAL GROUPS
+
+// Trailing edge and wake tip
+Physical Line("wakeTip") = {162};
+Physical Line("te") = {61};
+
+// Internal Field:
+Physical Volume("field") = {1};
+Physical Volume("field") += {2};
+
+// Wing:
+Physical Surface("wing") = {11,12,13,71,72,73};
+Physical Surface("wing_") = {14,15,16,74,75,76};
+Physical Surface("wing_upper") = {11,12,13,71,72,73};
+Physical Surface("wing_lower") = {14,15,16,74,75,76};
+
+// Tip
+//Physical Surface("tip") = {71,72,73,74,75,76};
+
+// Symmetry:
+Physical Surface("symmetry") = {1};
+Physical Surface("symmetry") += {6};
+
+// Farfield:
+Physical Surface("upstream") = {10};
+Physical Surface("farfield") = {3,4,5,8,9};
+
+// Downstream:
+Physical Surface("downstream") = {2};
+Physical Surface("downstream") += {7};
+
+// Wake:
+Physical Surface("wake") = {81};
+
+Coherence;
diff --git a/blast/tests/dart/onera_3D.py b/blast/tests/dart/onera_3D.py
new file mode 100644
index 0000000000000000000000000000000000000000..2ad42d166fd0c67413dcff7a633233836ee4d034
--- /dev/null
+++ b/blast/tests/dart/onera_3D.py
@@ -0,0 +1,157 @@
+#!/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()