merge Amaury.Bilocoq/waves (amaury) into feature_vii + cleaning. Tests to be fixed.

Squashed commit of the following:

commit c6c791de
Author: AmauryBilocq <amaury.bilocq@gmail.com>
Date:   Fri Nov 6 12:03:29 2020 +0100

    Quasi-simultaneous coupler for VII with test cases

commit 0cbee7e4
Merge: fd99999c ca62dcef
Author: AmauryBilocq <amaury.bilocq@gmail.com>
Date:   Sun Apr 26 14:58:26 2020 +0200

    V1.0 Viscous Inviscid interaction with direct solver

commit fd99999c
Author: AmauryBilocq <amaury.bilocq@gmail.com>
Date:   Sun Apr 26 14:47:01 2020 +0200

    V1.0 Viscous Inviscid interaction with direct solver

commit ca62dcef
Author: AmauryBilocq <amaury.bilocq@gmail.com>
Date:   Wed Oct 30 18:41:17 2019 +0100

    Trial to find the stagnation point

flow/default.py

@@ -65,11 +65,14 @@ def problem(msh, dim, alpha, beta, minf, mcrit, sref, lref, xref, yref, zref, su
     bnd = flo.Boundary(msh, [sur, fld])
     pbl.add(bnd)
     if vsc:
-        blw = flo.Blowing(msh, 'airfoil')
+        blw = flo.Blowing(msh, sur)
+        blww = flo.Blowing(msh, wk)
         pbl.add(blw)
+        pbl.add(blww)
     else:
         blw = None
-    return pbl, dirichlet, wake, bnd, blw
+        blww = None
+    return pbl, dirichlet, wake, bnd, [blw, blww]
 
 ## Initialize Picard solver
 def picard(pbl):

flow/tests/bli.py

@@ -16,9 +16,18 @@
 # limitations under the License.
 
 
-# @brief Compute lifting (linear or nonlinear) viscous flow around a naca0012
-# @authors Adrien Crovato, Amaury Bilocq
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+#
+# Amaury Bilocq
 # Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012: 
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, p = 2, m = 7, tol = 10^-4, msTE = 0.01, msLE = 0.001
+#    -> nIt = 41, Cl = 0.58, Cd = 0.0062, xtrTop = 0.0555, xtrBot = 0.7397
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, p = 2, m = 7, tol = 10^-4, msTE = 0.01, msLE = 0.001
+#    -> nIt = , Cl = 0.69, Cd = 0.0067, xtrTop = 0.0384, xtrBot = 0.7364
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, p = 2, m = 11, tol = 5.2*10^-3, msTE = 0.01, msLE = 0.00075
+#    -> nIt = , Cl = 1.39 , Cd = 0.011, xtrTop = 0.008, xtrBot = 1
 #
 # CAUTION
 # This test is provided to ensure that the solver works properly.
@@ -29,10 +38,11 @@ import flow.utils as floU
 import flow.default as floD
 import flow.viscous.solver as floVS
 import flow.viscous.coupler as floVC
+import tbox
 import tbox.utils as tboxU
 import fwk
 from fwk.testing import *
-from fwk.coloring import ccolors 
+from fwk.coloring import ccolors
 
 def main():
     # timer
@@ -40,16 +50,23 @@ def main():
     tms['total'].start()
 
     # define flow variables
-    alpha = 0*math.pi/180
-    M_inf = 0.0
+    Re = 1e7
+    alpha = 5*math.pi/180
+    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+    M_inf = 0.5
     M_crit = 5 # Squared critical Mach number (above which density is modified)
     c_ref = 1
     dim = 2
 
+    # define filter parameters and tolerance of the VII 
+    p = 2
+    m = 7
+    tol = 1e-4
+
     # mesh the geometry
     print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
     tms['msh'].start()
-    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.0075, 'msLe' : 0.0075}
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
@@ -62,8 +79,8 @@ def main():
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
     tms['solver'].start()
     isolver = floD.newton(pbl)
-    vsolver = floVS.Solver(blw)
-    coupler = floVC.Coupler(isolver, vsolver, gmshWriter)
+    vsolver = floVS.Solver(Re, p, m)
+    coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
     coupler.run()
     tms['solver'].stop()
 
@@ -72,8 +89,8 @@ def main():
     tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
     # display results
     print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
-    print('       M    alpha       Cl       Cd       Cm')
-    print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alpha*180/math.pi, isolver.Cl, isolver.Cd, isolver.Cm))
+    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(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
 
     # display timers
     tms['total'].stop()
@@ -83,16 +100,32 @@ def main():
 
     # visualize solution and plot results
     floD.initViewer(pbl)
-    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, isolver.Cd, isolver.Cm, 4), True)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
 
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
-    if M_inf == 0 and alpha == 0*math.pi/180:
-        tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.41, 1e-1)) # TODO check value and tolerance
-        tests.add(CTest('Cl', isolver.Cl, 0., 5e-2))
+    if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 0.58, 5e-2))
+        tests.add(CTest('Cd', vsolver.Cd, 0.0062, 0.01))
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 0.01))
+        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.056, 0.05))
+        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.05))
+    elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 0.69, 5e-2))
+        tests.add(CTest('Cd', vsolver.Cd, 0.0067, 0.01))
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 0.01))
+        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.038, 0.05))
+        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.05))
+    elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 1.39, 5e-2))
+        tests.add(CTest('Cd', vsolver.Cd, 0.0011, 0.01))
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 0.01))
+        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.008, 0.05))
+        tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.05))
     else:
         raise Exception('Test not defined for this flow')
+    
     tests.run()
 
     # eof

flow/viscous/airfoil.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 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.
+
+
+## Airfoil around which the boundary layer is computed
+#
+# Amaury Bilocq
+
+from flow.viscous.boundary import Boundary
+import numpy as np 
+
+class Airfoil(Boundary): 
+    def __init__(self, _boundary):
+        Boundary.__init__(self, _boundary)
+        self.name = 'airfoil' # type of boundary
+        self.T0 = 0 # initial condition for the momentum thickness
+        self.H0 = 0 
+        self.n0 = 0 
+        self.Ct0 = 0
+        self.TEnd = [0,0]
+        self.HEnd = [0,0]
+        self.CtEnd = [0,0]
+        self.nEnd = [0,0]
+
+    def initialConditions(self, Re, dv):
+        if dv > 0:
+            self.T0 = np.sqrt(0.075/(Re*dv))
+        else:
+            self.T0 = 1e-6
+        self.H0 = 2.23 # One parameter family Falkner Skan
+        self.n0 = 0
+        self.Ct0 = 0
+        return self.T0, self.H0, self.n0, self.Ct0
+
+    def connectList(self):
+        ''' Sort the value read by the viscous solver/ Create list of connectivity  
+        '''
+        N1 = np.zeros(self.nN, dtype=int) # Node number
+        connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
+        connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
+        data = np.zeros((self.boundary.nodes.size(), 10))
+        i = 0
+        for n in self.boundary.nodes:
+            data[i,0] = n.no
+            data[i,1] = n.pos[0]
+            data[i,2] = n.pos[1]
+            data[i,3] = n.pos[2]
+            data[i,4] = self.v[i,0]
+            data[i,5] = self.v[i,1]
+            data[i,6] = self.Me[i]
+            data[i,7] = self.rhoe[i]
+            data[i,8] = self.deltaStar[i]
+            data[i,9] = self.xx[i]
+            i += 1
+        # Table containing the element and its nodes
+        eData = np.zeros((self.nE,3), dtype=int) 
+        for i in range(0, self.nE):
+            eData[i,0] = self.boundary.tag.elems[i].no
+            eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
+            eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
+        # Defining TE/LE nodes number 
+        idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
+        globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
+        idxTE = np.where(data[:,1] == 1.0)[0]
+        if idxTE[0] < idxTE[1]:
+            upperIdxTE = idxTE[0]
+            lowerIdxTE = idxTE[1]
+        else:
+            upperIdxTE = idxTE[1]
+            lowerIdxTE = idxTE[0]
+        upperGlobTE = data[upperIdxTE,0] # Number of the upper TE node
+        lowerGlobTE = data[lowerIdxTE,0] # Number of the lower TE node   
+        connectListElems[0] = np.where(eData[:,1] == globStag)[0]
+        N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes
+        connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
+        i = 1
+        upperTE = 0 
+        lowerTE = 0
+        # Sort the suction part 
+        while upperTE == 0:
+            N1[i] = eData[connectListElems[i-1],2] # Second node of the element
+            connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes 
+            connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
+            if eData[connectListElems[i],2] == int(upperGlobTE): 
+                upperTE = 1
+            i += 1
+        # Sort the pressure side
+        connectListElems[i] = np.where(eData[:,2] == globStag)[0]
+        connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
+        N1[i] = eData[connectListElems[i],2]
+        while lowerTE == 0:
+            N1[i+1]  = eData[connectListElems[i],1] # First node of the element
+            connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes 
+            connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
+            if eData[connectListElems[i+1],1] == int(lowerGlobTE): 
+                lowerTE = 1
+            i += 1
+        connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
+        data[:,0] = data[connectListNodes,0]
+        data[:,1] = data[connectListNodes,1]
+        data[:,2] = data[connectListNodes,2]
+        data[:,3] = data[connectListNodes,3]
+        data[:,4] = data[connectListNodes,4]
+        data[:,5] = data[connectListNodes,5]
+        data[:,6] = data[connectListNodes,6]
+        data[:,7] = data[connectListNodes,7]
+        data[:,8] = data[connectListNodes,8]
+        data[:,9] = data[connectListNodes,9]
+        # Separated upper and lower part 
+        data = np.delete(data,0,1)       
+        uData = data[0:np.argmax(data[:,0])+1]
+        lData = data[np.argmax(data[:,0])+1:None]
+        lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
+        return connectListNodes, connectListElems,[uData, lData]

flow/viscous/boundary.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 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.
+
+
+## Base class representing a physical boundary
+#
+# Amaury Bilocq
+
+import numpy as np
+
+class Boundary: 
+    def __init__(self, _boundary):
+        self.boundary = _boundary # boundary of interest
+        self.nN = len(self.boundary.nodes) # number of nodes
+        self.nE = len(self.boundary.tag.elems) # number of elements
+        self.v = np.zeros((self.nN, 3)) # velocity at edge of BL
+        self.u = np.zeros(self.nE) # blowing velocity
+        self.Me = np.zeros(self.nN) # Mach number at the edge of the boundary layer 
+        self.rhoe = np.zeros(self.nN) # density at the edge of the boundary layer
+        self.connectListnodes = np.zeros(self.nN) # connectivity for nodes
+        self.connectListElems = np.zeros(self.nE) # connectivity for elements
+        self.deltaStar = np.zeros(self.nN) # displacement thickness
+        self.xx = np.zeros(self.nN) # coordinates in the reference of the physical surface. Used to store the values for the interaction law

flow/viscous/coupler.py

@@ -16,18 +16,21 @@
 # limitations under the License.
 
 
-# @brief Viscous-inviscid coupler
-# @authors Adrien Crovato, Amaury Bilocq
+## Viscous-inviscid coupler (quasi-simultaneous coupling)
+#
+# Amaury Bilocq
 
 import numpy as np
 from fwk.coloring import ccolors
+from flow.viscous.airfoil import Airfoil
+from flow.viscous.wake import Wake
 
 class Coupler:
-    def __init__(self, _isolver, _vsolver, _writer):
-        '''...
-        '''
+    def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer):
         self.isolver =_isolver # inviscid solver
         self.vsolver = _vsolver # viscous solver
+        self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects
+        self.tol = _tol # tolerance of the VII
         self.writer = _writer
 
     def run(self):
@@ -35,28 +38,48 @@ class Coupler:
         '''
         # initialize loop
         it = 0
-        converged = False # temp
-        while True:
+        converged = 0 # temp
+        CdOld = self.vsolver.Cd
+        while converged == 0:
             print(ccolors.ANSI_BLUE + 'Iteration: ', it, ccolors.ANSI_RESET)
             # run inviscid solver
             self.isolver.run()
-            # get velocity at edge of BL from inviscid solver and update viscous solver
-            for i in range(0, len(self.vsolver.boundary.nodes)):
-                self.vsolver.v[i,0] = self.isolver.U[self.vsolver.boundary.nodes[i].row][0]
-                self.vsolver.v[i,1] = self.isolver.U[self.vsolver.boundary.nodes[i].row][1]
-                self.vsolver.v[i,2] = self.isolver.U[self.vsolver.boundary.nodes[i].row][2]
-            # run viscous solver
-            self.vsolver.run()
-            # get blowing velocity from viscous solver and update inviscid solver
-            for i in range(0, self.vsolver.nE):
-                self.vsolver.boundary.setU(i, self.vsolver.u[i])
-            # convergence test
-            if converged:
-                break
+            print('--- Viscous solver parameters ---')
+            print('Filter windows length:', self.vsolver.m)
+            print('Filter polynomial order:', self.vsolver.p)
+            print('Tolerance:', self.tol)
+            print('--- Viscous problem definition ---')
+            print('Reynolds number:', self.vsolver.Re)
+            print('')
+            for n in range(0, len(self.group)):
+                print('Computing for', self.group[n].name, '...', end = ' ')
+                for i in range(0, len(self.group[n].boundary.nodes)):
+                    self.group[n].v[i,0] = self.isolver.U[self.group[n].boundary.nodes[i].row][0]
+                    self.group[n].v[i,1] = self.isolver.U[self.group[n].boundary.nodes[i].row][1]
+                    self.group[n].v[i,2] = self.isolver.U[self.group[n].boundary.nodes[i].row][2]
+                    self.group[n].Me[i] = self.isolver.M[self.group[n].boundary.nodes[i].row]
+                    self.group[n].rhoe[i] = self.isolver.rho[self.group[n].boundary.nodes[i].row]
+                    # run viscous solver
+                self.vsolver.run(self.group[n])
+                if self.group[n].name == 'airfoil':
+                    self.group[n+1].T0 = self.group[n].TEnd[0]+self.group[n].TEnd[1]
+                    self.group[n+1].H0 = (self.group[n].HEnd[0]*self.group[n].TEnd[0]+self.group[n].HEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0
+                    self.group[n+1].Ct0 = (self.group[n].CtEnd[0]*self.group[n].TEnd[0]+self.group[n].CtEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0
+                # get blowing velocity from viscous solver and update inviscid solver                    
+                for i in range(0, self.group[n].nE):                   
+                    self.group[n].boundary.setU(i, self.group[n].u[i])
+                print('done!')
+            print('    Iter           Cd      xtr_top      xtr_bot        Error')
+            print('{0:8d} {1:12.5f} {2:12.5f} {3:12.5f} {4:12.5f}'.format(it, self.vsolver.Cd, self.vsolver.xtr[0], self.vsolver.xtr[1], abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd))
+            print('')
+            # Converged or not
+            if abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd < self.tol:
+                converged = 1   
             else:
-                converged = True # perform 2 iterations
-                it += 1
-            
+                CdOld = self.vsolver.Cd 
+            it += 1
+            self.vsolver.it += 1     
         # save results
         self.isolver.save(0, self.writer)
+        self.vsolver.writeFile()
         print('\n')

flow/viscous/newton.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 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.
+
+
+## Newton raphson method coupled with linear solver
+#
+# Amaury Bilocq
+
+import numpy as np
+from fwk.coloring import ccolors
+
+def newton(f, x, maxIt, tol=1.0e-9):
+    '''Newton procedure to solve the non linear set of  boundary layer equations.
+    Boundary layer equations are parabolic in x -> downstream marching is used
+    An implicit marching model is applied (Crank Nicolson) -> matrix inversion is performed
+    '''
+    # Compute Jacobian 
+    def jacobian(f, x):
+        dx = 1.0e-8
+        n = len(x)
+        jac = np.zeros((n,n))
+        f0 = f(x) 
+        for j in range(n):
+            Dxj = (abs(x[j])*dx if x[j] != 0 else dx)
+            x_plus = [(xi if k != j else xi+Dxj) for k, xi in enumerate(x)]
+            jac[:,j] = (f(x_plus)-f0)/Dxj # numerical solution of the jacobian 
+        return jac, f0
+    # Solve
+    for i in range(maxIt):
+        jac, f0 = jacobian(f,x)
+        if np.sqrt(np.dot(f0,f0)/len(x)) < tol and x[0]>0 and  x[1]>0: return x 
+        dx = np.linalg.solve(jac, -f0)
+        x = x + dx
+        x[1] = max(x[1],1.0005) # to avoid too low shape parameter. Solver must be redone to be more robust
+        if np.sqrt(np.dot(dx, dx)) < tol*max(abs(any(x)),1) and x[0]>0 and x[1]>0: return x
+    # Error
+    raise RuntimeError(ccolors.ANSI_RED + 'Viscous Newton solver - too many iterations, aborting!\n' + ccolors.ANSI_RESET)

flow/viscous/wake.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 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.
+
+
+## Wake behind airfoil (around which the boundary layer is computed)
+#
+# Amaury Bilocq
+
+from flow.viscous.boundary import Boundary
+from flow.viscous.airfoil import Airfoil
+import numpy as np 
+
+class Wake(Boundary): 
+    def __init__(self, _boundary):
+        Boundary.__init__(self, _boundary)
+        self.name = 'wake' # type of boundary
+        self.T0 = 0 # initial condition for the momentum thickness
+        self.H0 = 0 
+        self.n0 = 9 # wake is always turbulent 
+        self.Ct0 = 0
+
+
+    def connectList(self):
+        ''' Sort the value read by the viscous solver/ Create list of connectivity  
+        '''
+        N1 = np.zeros(self.nN, dtype=int) # Node number
+        connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
+        connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
+        data = np.zeros((self.boundary.nodes.size(), 10))
+        i = 0
+        for n in self.boundary.nodes:
+            data[i,0] = n.no
+            data[i,1] = n.pos[0]
+            data[i,2] = n.pos[1]
+            data[i,3] = n.pos[2]
+            data[i,4] = self.v[i,0]
+            data[i,5] = self.v[i,1]
+            data[i,6] = self.Me[i]
+            data[i,7] = self.rhoe[i]
+            data[i,8] = self.deltaStar[i]
+            data[i,9] = self.xx[i]
+            i += 1
+        # Table containing the element and its nodes
+        eData = np.zeros((self.nE,3), dtype=int) 
+        for i in range(0, self.nE):
+            eData[i,0] = self.boundary.tag.elems[i].no
+            eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
+            eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
+        connectListNodes = data[:,1].argsort()
+        connectListElems[0] = np.where(eData[:,1] == 1.0)[0]
+        for i in range(1, len(eData[:,0])):
+            connectListElems[i] = np.where(eData[:,1] == eData[connectListElems[i-1],2])[0]
+        data[:,0] = data[connectListNodes,0]
+        data[:,1] = data[connectListNodes,1]
+        data[:,2] = data[connectListNodes,2]
+        data[:,3] = data[connectListNodes,3]
+        data[:,4] = data[connectListNodes,4]
+        data[:,5] = data[connectListNodes,5]
+        data[:,6] = data[connectListNodes,6]
+        data[:,7] = data[connectListNodes,7]
+        data[:,8] = data[connectListNodes,8]
+        data[:,9] = data[connectListNodes,9]
+        # Separated upper and lower part 
+        data = np.delete(data,0,1)
+        return connectListNodes, connectListElems, data