(feat) Sections can be generated using vtk

blast/blCoupler.py

@@ -53,6 +53,8 @@ class Coupler:
         print(f'{"Number of bodies:":<23s} {self.isol.getnBodies():<5.0f}')
         print(f'{"Number of sections:":<23s} {len(self.isol.cfg["sections"][0]):<5.0f}')
         print(f'{"Interpolator:":<23s} {self.isol.cfg["interpolator"]}')
+        if 'genSections' in self.isol.cfg:
+            print(f'{"Section generation:":<23s} {self.isol.cfg["genSections"]}')
         print('')
         print('-- Coupler parameters --') 
         print(f'{"Maximum number of iterations:":<30s} {self.maxIter:<3.0f}')
@@ -62,7 +64,6 @@ class Coupler:
         print(f'{"Verbosity level:":<30s} {self.isol.getVerbose():<2.0f}')
         print('')
 
-
     def run(self, write=True):
         # Aerodynamic coefficients.
         aeroCoeffs = {'Cl':[], 'Cd': [], 'Cdwake': []}

blast/interfaces/blSolversInterface.py

@@ -72,19 +72,25 @@ class SolversInterface:
             from blast.interfaces.interpolators.blRbfInterpolator import RbfInterpolator as interp
             print('Initializing RBF interpolator.')
             # Viscous sections
-            if 'vMsh' in self.cfg:
+            if self.cfg['genSections'] == 'file':
+                print('Getting sections from file')
+                if 'vMsh' not in self.cfg:
+                    raise RuntimeError('No viscous mesh provided.')
                 # If the mesh is provided we use it
                 print('Using provided viscous mesh')
                 self.getVWing()
-            else:
-                # Else generate the sections (upper and lower
-                # sides must be identified)
-                print('Generating sections')
+            elif self.cfg['genSections'] == 'auto':
+                print('Generating sections in auto mode')
                 for ibody in range(self.getnBodies()):
                     upNodes, lwNodes, upElms, lwElms = self.getUpperLowerSides()
                     self.iBnd[ibody][0].setUpperLower(upNodes[ibody], lwNodes[ibody], upElms[ibody], lwElms[ibody])
                 self.cfg['sidesIdentified'] = True
                 self.createVWing()
+            elif self.cfg['genSections'] == 'vtk':
+                print('Generating sections in vtk mode')
+                self.createVWingVTK()
+            else:
+                raise RuntimeError(f'Incorrect section generation method {self.cfg["genSections"]}')
         else:
             raise RuntimeError('Incorrect interpolator specified in config.')
 
@@ -307,7 +313,7 @@ class SolversInterface:
         ndim = self.getnDim()
         sections = self.cfg['sections']
         # Cosine spacing distribution
-        zeta = np.linspace(0., np.pi, self.cfg['nPoints'])
+        zeta = np.linspace(0., np.pi, self.cfg['nPoints']//2+1)
 
         for ibody in range(self.getnBodies()):
 
@@ -326,11 +332,11 @@ class SolversInterface:
             interpolated_upper = []
             interpolated_lower = []
             for i, y in enumerate(sections[ibody]):
-                data_up = np.zeros((self.cfg['nPoints'], 3))
-                data_lw = np.zeros((self.cfg['nPoints'], 3))
+                data_up = np.zeros((self.cfg['nPoints']//2+1, 3))
+                data_lw = np.zeros((self.cfg['nPoints']//2+1, 3))
                 # Select only few neighbors on inviscid mesh (10% of the span is kept)
-                local_upper_nodes = upper_nodes[abs(upper_nodes[:,1] - y) < .07 * (np.max(upper_nodes[:,1]) - np.min(upper_nodes[:,1])), :]
-                local_lower_nodes = lower_nodes[abs(lower_nodes[:,1] - y) < .07 * (np.max(lower_nodes[:,1]) - np.min(lower_nodes[:,1])), :]
+                local_upper_nodes = upper_nodes[abs(upper_nodes[:,1] - y) < .12 * (np.max(upper_nodes[:,1]) - np.min(upper_nodes[:,1])), :]
+                local_lower_nodes = lower_nodes[abs(lower_nodes[:,1] - y) < .12 * (np.max(lower_nodes[:,1]) - np.min(lower_nodes[:,1])), :]
 
                 # Compute the local chord to add extra points on viscous mesh
                 local_xmin_upper = np.min(local_upper_nodes[:,0])
@@ -360,7 +366,7 @@ class SolversInterface:
                     dummy[-1,ndim-1] = dummy[0, ndim-1]
                     # Smooth leading edge
                     mask = dummy[:,0] < (0.1*(np.max(dummy[:,0]) - np.min(dummy[:,0])) + np.min(dummy[:,0]))
-                    spl_n, u_n = make_splprep([dummy[mask,0], dummy[mask,2]], s=1e-8)
+                    spl_n, u_n = make_splprep([dummy[mask,0], dummy[mask,2]], s=1e-7)
                     dummy[mask, 0] = spl_n(u_n)[0]
                     dummy[mask, 2] = spl_n(u_n)[1]
                     data_up = dummy[:data_up.shape[0], :]
@@ -370,13 +376,13 @@ class SolversInterface:
             # Create sections
             # for i in range(len(interpolated_upper)):
                 chord_upper = np.max(interpolated_upper[i][:,0]) - np.min(interpolated_upper[i][:,0])
-                data_up = np.zeros((self.cfg['nPoints'], 3))
+                data_up = np.zeros((self.cfg['nPoints']//2+1, 3))
                 data_up[:,0] = ( chord_upper / 2) * (np.cos(zeta) + 1) + np.min(interpolated_upper[i][:,0])
                 data_up[:,1] = sections[ibody][i] * np.ones(data_up.shape[0])
                 data_up[:,ndim-1] = interp1d(interpolated_upper[i][:,0], interpolated_upper[i][:,ndim-1], kind='linear', fill_value='extrapolate')(data_up[:,0])
 
                 chord_lower = np.max(interpolated_lower[i][:,0]) - np.min(interpolated_lower[i][:,0])
-                data_lw = np.zeros((self.cfg['nPoints'], 3))
+                data_lw = np.zeros((self.cfg['nPoints']//2+1, 3))
                 data_lw[:,0] = (chord_lower / 2) * (np.cos(zeta) + 1) + np.min(interpolated_lower[i][:,0])
                 data_lw[:,1] = sections[ibody][i] * np.ones(data_lw.shape[0])
                 data_lw[:,ndim-1] = interp1d(interpolated_lower[i][:,0], interpolated_lower[i][:,ndim-1], kind='linear', fill_value='extrapolate')(data_lw[:,0])
@@ -438,6 +444,80 @@ class SolversInterface:
                 self.vBnd[ibody][iy][1].elemsCoord = elems_coords_wake[iy]
             print('Created', len(sections_final), 'sections and', len(wake_final), 'wake sections')
 
+    def createVWingVTK(self):
+        try:
+            import vtk
+        except:
+            raise Exception('VTK not found!\n')
+        import tboxVtk.reader as vtkR
+        import tboxVtk.cutter as vtkC
+        import tbox.utils as tU
+        reader = vtkR.Reader()
+
+        print('oui')
+        self.save(sfx="")
+        print('saved')
+
+        mshName = self.getMeshName()
+        reader.open(mshName)
+        cutter = vtkC.Cutter(reader.reader.GetOutputPort())
+
+        zeta = np.linspace(0., np.pi, self.cfg['nPoints']//2+1)
+
+        for ibody, ys in enumerate(self.cfg['sections']):
+            wingTag = self.getWingTag(ibody)
+            wakeTag = self.getWakeTag(ibody)
+            print('tag')
+            for isec in range(0, len(ys)):
+                pts, elems, _ = cutter.extract(cutter.cut(wingTag, [0., ys[isec], 0.], [0., 1., 0.]), 2, ['Cp'])
+                tU.sort(elems, pts)
+                pts = np.vstack((pts, pts[0,:])) # Duplicate TE point
+                ile = np.argmin(pts[:,0]) # LE index
+                ite = np.argmax(pts[:,0]) # TE index
+                upper_side = pts[ile:, :]
+                lower_side = pts[:ile+1, :]
+                # Remove duplicate x
+                _, unique_upper = np.unique(upper_side[:, 0], return_index=True)
+                unique_upper = np.sort(unique_upper)
+                upper_side = upper_side[unique_upper]
+                _, unique_lower = np.unique(lower_side[:, 0], return_index=True)
+                unique_lower = np.sort(unique_lower)
+                lower_side = lower_side[unique_lower]
+
+                chord = pts[ite, 0] - pts[ile, 0]
+                x_interp = (chord/2) * (np.cos(zeta)+1) + pts[ile,0]
+                upper_interp = interp1d(upper_side[:,0], upper_side[:,2], kind='linear')(x_interp)
+                lower_interp = interp1d(lower_side[:,0], lower_side[:,2], kind='linear')(x_interp)
+                upper_final = np.column_stack((x_interp, np.full_like(x_interp, ys[isec]), upper_interp))
+                lower_final = np.column_stack((x_interp, np.full_like(x_interp, ys[isec]), lower_interp))
+                section = np.row_stack((upper_final, np.flip(lower_final, axis=0)[1:]))
+
+                # Compute elements positions
+                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
+                
+                # Wake section
+                ptsw, _, _ = cutter.extract(cutter.cut(wakeTag, [0., ys[isec], 0.], [0., 1., 0.]), 2, ['Cp'])
+                # Sort in ascending order of the first columns
+                ptsw = ptsw[ptsw[:,0].argsort()]
+
+                chordw = ptsw[-1, 0] - ptsw[0, 0]
+                x_interp = np.flip((chordw/2) * (np.cos(zeta)+1) + ptsw[0,0])
+                interp = interp1d(ptsw[:,0], ptsw[:,2], kind='linear')(x_interp)
+                wake_final = np.column_stack((x_interp, np.full_like(x_interp, ys[isec]), interp))
+                
+                elems_coords_wake = np.empty((wake_final.shape[0]-1, wake_final.shape[1]), dtype=float, order='C')
+                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][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][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
+
     def getVWing(self):
         """Obtain the nodes' location and row and cg of all elements.
         If the mesh is the viscous one, TE nodes are duplicated (by creating a new row).

blast/interfaces/dart/blDartInterface.py

@@ -200,6 +200,15 @@ class DartInterface(SolversInterface):
 
     def getnBlowing(self, ibody=0):
         return len(self.blw[ibody])
+    
+    def getMeshName(self):
+        return self.iobj['msh'].name
+    
+    def getWingTag(self, ibody=0):
+        return self.blw[ibody][0].tag.no
+
+    def getWakeTag(self, ibody=0):
+        return self.blw[ibody][1].tag.no
 
     def reset(self):
         """Reset the solution.

blast/interfaces/interpolators/blRbfInterpolator.py

@@ -54,7 +54,7 @@ class RbfInterpolator(Interpolator):
                 Mi_wk = iDict[ibody][1].M
                 rhoi_wk = iDict[ibody][1].Rho
                 vi_wk = iDict[ibody][1].V[:,:self.ndim]
-                self.tms['getIdict'].start()
+                self.tms['getIdict'].stop()
 
                 self.tms['getVdict'].start()
                 # Viscous side
@@ -125,6 +125,7 @@ class RbfInterpolator(Interpolator):
                     sec[0].updateVars(v_airfoil, M_airfoil, rho_airfoil)
                     sec[1].updateVars(v_wake, M_wake, rho_wake)
                 self.tms['setting'].stop()
+            self.tms['tot'].stop()
         else:
             for ibody in range(len(iDict)):
                 for iSec in range(len(vDict[ibody])):
@@ -137,7 +138,6 @@ class RbfInterpolator(Interpolator):
                         M = self.__rbfToSection(iDict[ibody][iReg].nodesCoord[:,:(self.ndim if iDict[ibody][iReg].name == 'iWing' else 1)], iDict[ibody][iReg].M, reg.nodesCoord[:,:(self.ndim if 'vAirfoil' in reg.name else 1)])
                         rho = self.__rbfToSection(iDict[ibody][iReg].nodesCoord[:,:(self.ndim if iDict[ibody][iReg].name == 'iWing' else 1)], iDict[ibody][iReg].Rho, reg.nodesCoord[:,:(self.ndim if 'vAirfoil' in reg.name else 1)])
                         vDict[ibody][iSec][iReg].updateVars(v, M, rho)
-            self.tms['tot'].stop()
 
     def viscousToInviscid(self, iDict, vDict):
         

blast/tests/dart/t_rae2822_3D.py

@@ -78,14 +78,23 @@ def cfgInviscid(nthrds, verb):
 
 def cfgBlast(verb):
     return {
+        # General parameters
+        'writeSections': [0.2, 0.4, 0.6, 0.8, 1.0], # Spanwise locations to write the solution
+        'Sym':[0.],                                 # Symmetry plane
+        'spans': [1.],                              # Span of the wing
+
+        # Viscous flow parameters
         'Re' : 1e7,       # Freestream Reynolds number
         'Minf' : 0.8,     # Freestream Mach number (used for the computation of the time step only)
         'CFL0' : 1,       # Inital CFL number of the calculation
+        'xtrF' : [0., 0.],  # Forced transition locations
 
+        # Sections generation parameters
         'sections' : [np.linspace(0.05, 0.9, 10)],  # Sections on the wing
-        'writeSections': [0.2, 0.4, 0.6, 0.8, 1.0], # Spanwise locations to write the solution
-        'Sym':[0.],                                 # Symmetry plane
-        'spans': [1.],                              # Span of the wing
+        'genSections': 'auto',                       # Section generation technique ['auto', 'file', 'vtk']
+        'nPoints': 300,                             # Section will have nPoints+1 points because of the duplicate TE
+
+        # Interpolator parameters
         'interpolator': 'Rbf',                      # Interpolator type
         'rbftype': 'linear',                        # Radial basis function (rbf) type
         'smoothing': 1e-8,                          # rbf smoothing factor
@@ -93,13 +102,12 @@ def cfgBlast(verb):
         'neighbors': 10,                            # Number of neighbors for the interpolator
         'saveTag': 4,                               # Tag to save the solution with VTK
 
+        # Coupling parameters
         'Verb': verb,       # Verbosity level of the solver
-        'couplIter': 5,     # Maximum number of iterations
-        'couplTol' : 1e-2,  # Tolerance of the VII methodology
+        'couplIter': 50,     # Maximum number of iterations
+        'couplTol' : 1e-6,  # Tolerance of the VII methodology
         'iterPrint': 1,     # int, number of iterations between outputs
         'resetInv' : False, # bool, flag to reset the inviscid calculation at every iteration.
-        'xtrF' : [0., 0.],  # Forced transition locations
-        'nPoints': 300,
     }
 
 def main():