diff --git a/RAE2822_Airfoil_HISA/0_orig/T b/RAE2822_Airfoil_HISA/0_orig/T
new file mode 100644
index 0000000..ae3a909
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/0_orig/T
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#include         "include/freestreamConditions"
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform $T;
+
+boundaryField
+{
+    "wing.*"
+    {
+        type            characteristicWallTemperature;
+        value           $internalField;
+    }
+    symmetry1
+    {
+        type            symmetry;
+    }
+    symmetry2
+    {
+        type            symmetry;
+    }
+    inout
+    {
+        type            characteristicFarfieldTemperature;
+        #include        "include/freestreamConditions"
+        value           $internalField;
+        //type            inletOutlet;
+        //inletValue      $internalField;
+        //value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/0_orig/U b/RAE2822_Airfoil_HISA/0_orig/U
new file mode 100644
index 0000000..5122bab
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/0_orig/U
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include         "include/freestreamConditions"
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField uniform $U;
+
+boundaryField
+{
+    "wing.*"
+    {
+        type            boundaryCorrectedFixedValue;
+        value           uniform (0 0 0);
+    }
+    symmetry1
+    {
+        type            symmetry;
+    }
+    symmetry2
+    {
+        type            symmetry;
+    }
+    inout
+    {
+        type            characteristicFarfieldVelocity;
+        #include        "include/freestreamConditions"
+        value           uniform (0 0 0);
+        //type            inletOutlet;
+        //inletValue      $internalField;
+        //value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/0_orig/alphat b/RAE2822_Airfoil_HISA/0_orig/alphat
new file mode 100644
index 0000000..fdddca5
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/0_orig/alphat
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    "wing.*"
+    {
+        type            boundaryCorrectedFixedValue;
+        value           uniform 0;
+    }
+    symmetry1
+    {
+        type            symmetry;
+    }
+    symmetry2
+    {
+        type            symmetry;
+    }
+    inout
+    {
+        type            calculated;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/0_orig/include/freestreamConditions b/RAE2822_Airfoil_HISA/0_orig/include/freestreamConditions
new file mode 100644
index 0000000..98ce389
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/0_orig/include/freestreamConditions
@@ -0,0 +1,14 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+|                                                                             |
+|    HiSA: High Speed Aerodynamic solver                                      |
+|    Copyright (C) 2014-2017 Johan Heyns - CSIR, South Africa                 |
+|    Copyright (C) 2014-2017 Oliver Oxtoby - CSIR, South Africa               |
+|                                                                             |
+\*---------------------------------------------------------------------------*/
+
+U             (233.603 0 0);
+p             108988.0;
+T             255.56;
+#inputMode    merge
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/0_orig/nuTilda b/RAE2822_Airfoil_HISA/0_orig/nuTilda
new file mode 100644
index 0000000..b6e139b
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/0_orig/nuTilda
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nuTilda;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 4.5e-05;
+
+boundaryField
+{
+    "wing.*"
+    {
+        type            boundaryCorrectedFixedValue;
+        value           uniform 0.0;
+    }
+    symmetry1
+    {
+        type            symmetry;
+    }
+    symmetry2
+    {
+        type            symmetry;
+    }
+    inout
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/0_orig/nut b/RAE2822_Airfoil_HISA/0_orig/nut
new file mode 100644
index 0000000..9590a53
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/0_orig/nut
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    "wing.*"
+    {
+        type            boundaryCorrectedNutLowReWallFunction;
+        value           uniform 0;
+    }
+    symmetry1
+    {
+        type            symmetry;
+    }
+    symmetry2
+    {
+        type            symmetry;
+    }
+    inout
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+}
+
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/0_orig/p b/RAE2822_Airfoil_HISA/0_orig/p
new file mode 100644
index 0000000..0a2fb93
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/0_orig/p
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#include         "include/freestreamConditions"
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform $p;
+
+boundaryField
+{
+    "wing.*"
+    {
+        type            characteristicWallPressure;
+        value           $internalField;
+    }
+    symmetry1
+    {
+        type            symmetry;
+    }
+    symmetry2
+    {
+        type            symmetry;
+    }
+    inout
+    {
+        type            characteristicFarfieldPressure;
+        #include        "include/freestreamConditions"
+        value           $internalField;
+
+        //type            outletInlet;
+        //outletValue     $internalField;
+        //value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/Allclean.sh b/RAE2822_Airfoil_HISA/Allclean.sh
new file mode 100644
index 0000000..8bfb3cd
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/Allclean.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+
+while true 
+do
+    read -p "Delete everything and resume to the default setup (y/n)?" yn
+    case $yn in
+        [Yy]* ) 
+            # clean everyting
+            echo "Cleaning..."
+            rm -rf 0
+            rm -rf postProcessing
+            rm -rf constant/extendedFeatureEdgeMesh
+            rm -rf constant/triSurface
+            rm -rf constant/polyMesh/
+            rm -rf *.bin *.info *.dat *.xyz *.stl *.txt
+            rm -rf processor* 0.0000*
+            rm -rf {1..9}*
+            exit
+            ;;
+        [Nn]* ) exit;;
+        * ) echo "Please answer yes or no.";;
+    esac
+done
+
diff --git a/RAE2822_Airfoil_HISA/constant/thermophysicalProperties b/RAE2822_Airfoil_HISA/constant/thermophysicalProperties
new file mode 100644
index 0000000..c4a3e4c
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/constant/thermophysicalProperties
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+|                                                                             |
+|    HiSA: High Speed Aerodynamic solver                                      |
+|    Copyright (C) 2014-2017 Johan Heyns - CSIR, South Africa                 |
+|    Copyright (C) 2014-2017 Oliver Oxtoby - CSIR, South Africa               |
+|                                                                             |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            hePsiThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+
+mixture
+{
+   specie
+   {
+       nMoles          1;
+       molWeight       28.966;
+   }
+   thermodynamics
+   {
+       Cp              1005;
+       Hf              0;
+   }
+   transport
+   {
+        mu                  1.8e-5; 
+        Pr                  0.700000; 
+        TRef                300.000000; 
+   }
+}
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/constant/turbulenceProperties b/RAE2822_Airfoil_HISA/constant/turbulenceProperties
new file mode 100644
index 0000000..9aab251
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/constant/turbulenceProperties
@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+|                                                                             |
+|    HiSA: High Speed Aerodynamic solver                                      |
+|    Copyright (C) 2014-2017 Johan Heyns - CSIR, South Africa                 |
+|    Copyright (C) 2014-2017 Oliver Oxtoby - CSIR, South Africa               |
+|                                                                             |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RAS;
+
+RAS
+{
+    RASModel        SpalartAllmaras;
+
+    turbulence      on;
+
+    printCoeffs     on;
+Prt 0.85;
+}
+
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/genAirFoilMesh.py b/RAE2822_Airfoil_HISA/genAirFoilMesh.py
new file mode 100644
index 0000000..4b9bae1
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/genAirFoilMesh.py
@@ -0,0 +1,269 @@
+#!/usr/bin/env python
+"""
+This script reads a coarse airfoil profile, refine the profile using spline function,
+and outputs it as surfaceMesh.xyz. Then it generate a 3D volume mesh with nSpan layers
+in the z direction using pyHyp, available at https://github.com/mdolab/pyhyp
+Note: the airfoil data should be seperated into PS and SS surfaces, they should start from the
+LE and ends at TE. We use blunt TE so truncate the PS and SS data at about 99.8% of the chord.
+"""
+
+from pyhyp import pyHyp
+import numpy
+from pyspline import *
+
+########## user input ################
+# 2D
+prefix = "./profiles/"
+airfoilProfilePS = prefix + "RAE2822PS.profile"
+airfoilProfileSS = prefix + "RAE2822SS.profile"
+ZSpan = 0.01  # width in the z direction
+nSpan = 2  # how many points in z
+# PS parameters
+dX1PS = 0.0025  # first dx from the LE
+Alpha1PS = 1.2  # clustering from the LE
+dX2PS = 5e-4  # first dx from the TE
+Alpha2PS = 1.2  # clustering from the TE
+dXMaxPS = 0.015  # max dx for PS
+# SS parameters
+dX1SS = 0.0025  # first dx from the LE
+Alpha1SS = 1.2  # clustering from the LE
+dX2SS = 5e-4  # first dx from the TE
+Alpha2SS = 1.2  # clustering from the TE
+dXMaxSS = 0.015  # max dx for SS
+# TE parameters
+NpTE = 5  # number of points for blunt TE
+# 3D
+NpExtrude = 65  # how many points to extrude for the 3D volume mesh in y
+yWall = 1e-5  # first layer mesh length
+marchDist = 30.0  # march distance for extruding
+########## user input ################
+
+# read profiles
+fPS = open(airfoilProfilePS, "r")
+linesPS = fPS.readlines()
+fPS.close()
+xPS = []
+yPS = []
+zPS = []
+for line in linesPS:
+    cols = line.split()
+    xPS.append(float(cols[0]))
+    yPS.append(float(cols[1]))
+
+for i in range(len(xPS)):
+    zPS.append(0.0)
+
+fSS = open(airfoilProfileSS, "r")
+linesSS = fSS.readlines()
+fSS.close()
+xSS = []
+ySS = []
+zSS = []
+for line in linesSS:
+    cols = line.split()
+    xSS.append(float(cols[0]))
+    ySS.append(float(cols[1]))
+for i in range(len(xSS)):
+    zSS.append(0.0)
+
+# ------------PS
+# first compute how many stretching points do we need for both ends
+tmp = dX1PS
+for i in range(1000):
+    if tmp > dXMaxPS:
+        nStretch1PS = i
+        break
+    else:
+        tmp = tmp * Alpha1PS
+# print (nStretch1PS)
+tmp = dX2PS
+for i in range(1000):
+    if tmp > dXMaxPS:
+        nStretch2PS = i
+        break
+    else:
+        tmp = tmp * Alpha2PS
+# print (nStretch2PS)
+
+# now compute how much length does these two stretching end and the constant portion have
+xLPSConst = xPS[-1]
+xLPS1 = 0
+xLPS2 = 0
+for i in range(nStretch1PS):
+    xLPS1 += dX1PS * (Alpha1PS ** i)
+    xLPSConst -= dX1PS * (Alpha1PS ** i)
+for i in range(nStretch2PS):
+    xLPS2 += dX2PS * (Alpha2PS ** i)
+    xLPSConst -= dX2PS * (Alpha2PS ** i)
+# print xLPS1,xLPS2,xLPSConst
+# update dXMax, we just want to make sure these three portion add up
+nXConstPS = int(xLPSConst / dXMaxPS)
+dXMaxPS = xLPSConst / nXConstPS
+# print(dXMaxPS)
+
+# now we can add these three portions together
+xInterpPS = [0]
+tmp = dX1PS
+for i in range(nStretch1PS):
+    xInterpPS.append(xInterpPS[-1] + tmp)
+    tmp = tmp * Alpha1PS
+for i in range(nXConstPS):
+    xInterpPS.append(xInterpPS[-1] + dXMaxPS)
+tmp = dX2PS * (Alpha2PS ** (nStretch2PS - 1))
+for i in range(nStretch2PS):
+    xInterpPS.append(xInterpPS[-1] + tmp)
+    tmp /= Alpha2PS
+# print xInterpPS
+# Finally, we interpolate the refined stretch stuff
+c1PS = Curve(x=xPS, y=yPS, z=zPS, k=3)
+XPS = c1PS(xInterpPS)
+c2PS = Curve(X=XPS, k=3)
+x1PS = c2PS.X[:, 0]
+y1PS = c2PS.X[:, 1]
+
+
+# ------------SS
+# first compute how many stretching points do we need for both ends
+tmp = dX1SS
+for i in range(1000):
+    if tmp > dXMaxSS:
+        nStretch1SS = i
+        break
+    else:
+        tmp = tmp * Alpha1SS
+# print (nStretch1SS)
+tmp = dX2SS
+for i in range(1000):
+    if tmp > dXMaxSS:
+        nStretch2SS = i
+        break
+    else:
+        tmp = tmp * Alpha2SS
+# print (nStretch2SS)
+
+# now compute how much length does these two stretching end and the constant portion have
+xLSSConst = xSS[-1]
+xLSS1 = 0
+xLSS2 = 0
+for i in range(nStretch1SS):
+    xLSS1 += dX1SS * (Alpha1SS ** i)
+    xLSSConst -= dX1SS * (Alpha1SS ** i)
+for i in range(nStretch2SS):
+    xLSS2 += dX2SS * (Alpha2SS ** i)
+    xLSSConst -= dX2SS * (Alpha2SS ** i)
+# print xLSS1,xLSS2,xLSSConst
+# update dXMax, we just want to make sure these three portion add up
+nXConstSS = int(xLSSConst / dXMaxSS)
+dXMaxSS = xLSSConst / nXConstSS
+# print(dXMaxSS)
+
+# now we can add these three portions together
+xInterpSS = [0]
+tmp = dX1SS
+for i in range(nStretch1SS):
+    xInterpSS.append(xInterpSS[-1] + tmp)
+    tmp = tmp * Alpha1SS
+for i in range(nXConstSS):
+    xInterpSS.append(xInterpSS[-1] + dXMaxSS)
+tmp = dX2SS * (Alpha2SS ** (nStretch2SS - 1))
+for i in range(nStretch2SS):
+    xInterpSS.append(xInterpSS[-1] + tmp)
+    tmp /= Alpha2SS
+# print xInterpSS
+# Finally, we interpolate the refined stretch stuff
+c1SS = Curve(x=xSS, y=ySS, z=zSS, k=3)
+XSS = c1SS(xInterpSS)
+c2SS = Curve(X=XSS, k=3)
+x1SS = c2SS.X[:, 0]
+y1SS = c2SS.X[:, 1]
+
+# Since the TE is open we need to close it. Close it multiple linear segments.
+delta_y = numpy.linspace(y1PS[-1], y1SS[-1], NpTE, "d")
+delta_y = delta_y[1:]
+delta_x = numpy.ones_like(delta_y, "d")
+for i in range(len(delta_x)):
+    delta_x[i] = x1SS[-1]
+
+x1SS_Flip = x1SS[::-1]  # reverse the array #numpy.flip(x1SS,axis=0)
+xAll = numpy.append(x1SS_Flip, x1PS[1:])
+xAll = numpy.append(xAll, delta_x)
+
+y1SS_Flip = y1SS[::-1]  # reverse the array # numpy.flip(y1SS,axis=0)
+yAll = numpy.append(y1SS_Flip, y1PS[1:])
+yAll = numpy.append(yAll, delta_y)
+
+
+# print mesh statistics
+print("nPoints for PS: ", nStretch1PS + nStretch2PS + nXConstPS)
+print("nPoints for SS: ", nStretch1SS + nStretch2SS + nXConstSS)
+print("nPoints for TE: ", NpTE)
+print("nPoints Total: ", nStretch1PS + nStretch2PS + nXConstPS + nStretch1SS + nStretch2SS + nXConstSS + NpTE)
+print(
+    "Mesh cells: ",
+    (nStretch1PS + nStretch2PS + nXConstPS + nStretch1SS + nStretch2SS + nXConstSS + NpTE - 1)
+    * (NpExtrude - 1)
+    * (nSpan - 1),
+)
+
+# plt.plot(xPS,yPS,'-k',linewidth=1)
+# plt.plot(xAll,yAll,'ro',markersize=2)
+# plt.plot(xSS,ySS,'-k',linewidth=1)
+# plt.gca().set_aspect('equal', adjustable='box')
+# plt.show()
+# plt.savefig('figure.png',bbox_inches='tight')   # save the figure to file
+# plt.close()    # close the figure
+
+
+# Write the plot3d input file:
+f = open("surfaceMesh.xyz", "w")
+f.write("1\n")
+f.write("%d %d %d\n" % (len(xAll), nSpan, 1))
+for iDim in range(3):
+    for z in numpy.linspace(0.0, ZSpan, nSpan):
+        for i in range(len(xAll)):
+            if iDim == 0:
+                f.write("%20.16f\n" % xAll[i])
+            elif iDim == 1:
+                f.write("%20.16f\n" % yAll[i])
+            else:
+                f.write("%20.16f\n" % z)
+f.close()
+
+options = {
+    # ---------------------------
+    #        Input Parameters
+    # ---------------------------
+    "inputFile": "surfaceMesh.xyz",
+    "unattachedEdgesAreSymmetry": False,
+    "outerFaceBC": "farfield",
+    "autoConnect": True,
+    "BC": {1: {"jLow": "zSymm", "jHigh": "zSymm"}},
+    "families": "wall",
+    # ---------------------------
+    #        Grid Parameters
+    # ---------------------------
+    "N": NpExtrude,
+    "s0": yWall,
+    "marchDist": marchDist,
+    # ---------------------------
+    #   Pseudo Grid Parameters
+    # ---------------------------
+    "ps0": -1.0,
+    "pGridRatio": -1.0,
+    "cMax": 1.0,
+    # ---------------------------
+    #   Smoothing parameters
+    # ---------------------------
+    "epsE": 2.0,
+    "epsI": 4.0,
+    "theta": 2.0,
+    "volCoef": 0.20,
+    "volBlend": 0.0005,
+    "volSmoothIter": 20,
+}
+
+
+hyp = pyHyp(options=options)
+hyp.run()
+# hyp.writeCGNS('volumeMesh.cgns')
+hyp.writePlot3D("volumeMesh.xyz")
diff --git a/RAE2822_Airfoil_HISA/paraview.foam b/RAE2822_Airfoil_HISA/paraview.foam
new file mode 100644
index 0000000..e69de29
diff --git a/RAE2822_Airfoil_HISA/preProcessing.sh b/RAE2822_Airfoil_HISA/preProcessing.sh
new file mode 100644
index 0000000..2d4744c
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/preProcessing.sh
@@ -0,0 +1,20 @@
+#!/bin/bash
+
+if [ -z "$WM_PROJECT" ]; then
+  echo "OpenFOAM environment not found, forgot to source the OpenFOAM bashrc?"
+  exit
+fi
+
+# pre-processing
+
+# generate mesh
+echo "Generating mesh.."
+python genAirFoilMesh.py &> logMeshGeneration.txt
+plot3dToFoam -noBlank volumeMesh.xyz >> logMeshGeneration.txt
+autoPatch 45 -overwrite >> logMeshGeneration.txt
+createPatch -overwrite >> logMeshGeneration.txt
+renumberMesh -overwrite >> logMeshGeneration.txt
+echo "Generating mesh.. Done!"
+
+# copy initial and boundary condition files
+cp -r 0_orig 0
diff --git a/RAE2822_Airfoil_HISA/profiles/RAE2822PS.profile b/RAE2822_Airfoil_HISA/profiles/RAE2822PS.profile
new file mode 100644
index 0000000..bddfd79
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/profiles/RAE2822PS.profile
@@ -0,0 +1,63 @@
+  0.000000  0.000000
+  0.000602  -.003160
+  0.002408  -.006308
+  0.005412  -.009443
+  0.009607  -.012559
+  0.014984  -.015649
+  0.021530  -.018707
+  0.029228  -.021722
+  0.038060  -.024685
+  0.048005  -.027586
+  0.059039  -.030416
+  0.071136  -.033170
+  0.084265  -.035843
+  0.098396  -.038431
+  0.113495  -.040929
+  0.129524  -.043326
+  0.146447  -.045610
+  0.164221  -.047773
+  0.182803  -.049805
+  0.202150  -.051694
+  0.222215  -.053427
+  0.242949  -.054994
+  0.264302  -.056376
+  0.286222  -.057547
+  0.308658  -.058459
+  0.331555  -.059046
+  0.354858  -.059236
+  0.378510  -.058974
+  0.402455  -.058224
+  0.426635  -.056979
+  0.450991  -.055257
+  0.475466  -.053099
+  0.500000  -.050563
+  0.524534  -.047719
+  0.549009  -.044642
+  0.573365  -.041397
+  0.597545  -.038043
+  0.621490  -.034631
+  0.645142  -.031207
+  0.668445  -.027814
+  0.691342  -.024495
+  0.713778  -.021289
+  0.735698  -.018232
+  0.757051  -.015357
+  0.777785  -.012690
+  0.797850  -.010244
+  0.817197  -.008027
+  0.835779  -.006048
+  0.853553  -.004314
+  0.870476  -.002829
+  0.886505  -.001592
+  0.901604  -.000600
+  0.915735  0.000157
+  0.928864  0.000694
+  0.940961  0.001033
+  0.951995  0.001197
+  0.961940  0.001212
+  0.970772  0.001112
+  0.978470  0.000935
+  0.985016  0.000719
+  0.990393  0.000497
+  0.994588  0.000296
+  0.997592  0.000137
\ No newline at end of file
diff --git a/RAE2822_Airfoil_HISA/profiles/RAE2822SS.profile b/RAE2822_Airfoil_HISA/profiles/RAE2822SS.profile
new file mode 100644
index 0000000..bc472c9
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/profiles/RAE2822SS.profile
@@ -0,0 +1,63 @@
+  0.000000  0.000000
+  0.000602  0.003165
+  0.002408  0.006306
+  0.005412  0.009416
+  0.009607  0.012480
+  0.014984  0.015489
+  0.021530  0.018441
+  0.029228  0.021348
+  0.038060  0.024219
+  0.048005  0.027062
+  0.059039  0.029874
+  0.071136  0.032644
+  0.084265  0.035360
+  0.098396  0.038011
+  0.113495  0.040585
+  0.129524  0.043071
+  0.146447  0.045457
+  0.164221  0.047729
+  0.182803  0.049874
+  0.202150  0.051885
+  0.222215  0.053753
+  0.242949  0.055470
+  0.264302  0.057026
+  0.286222  0.058414
+  0.308658  0.059629
+  0.331555  0.060660
+  0.354858  0.061497
+  0.378510  0.062133
+  0.402455  0.062562
+  0.426635  0.062779
+  0.450991  0.062774
+  0.475466  0.062530
+  0.500000  0.062029
+  0.524534  0.061254
+  0.549009  0.060194
+  0.573365  0.058845
+  0.597545  0.057218
+  0.621490  0.055344
+  0.645142  0.053258
+  0.668445  0.050993
+  0.691342  0.048575
+  0.713778  0.046029
+  0.735698  0.043377
+  0.757051  0.040641
+  0.777785  0.037847
+  0.797850  0.035017
+  0.817197  0.032176
+  0.835779  0.029347
+  0.853553  0.026554
+  0.870476  0.023817
+  0.886505  0.021153
+  0.901604  0.018580
+  0.915735  0.016113
+  0.928864  0.013769
+  0.940961  0.011562
+  0.951995  0.009508
+  0.961940  0.007622
+  0.970772  0.005915
+  0.978470  0.004401
+  0.985016  0.003092
+  0.990393  0.002001
+  0.994588  0.001137
+  0.997592  0.000510
\ No newline at end of file
diff --git a/RAE2822_Airfoil_HISA/runScript.py b/RAE2822_Airfoil_HISA/runScript.py
new file mode 100644
index 0000000..bbebf99
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/runScript.py
@@ -0,0 +1,98 @@
+#!/usr/bin/env python
+"""
+DAFoam run script for the RAE2822 airfoil at transonic conditions
+"""
+
+# =============================================================================
+# Imports
+# =============================================================================
+import os
+import argparse
+import numpy as np
+from mpi4py import MPI
+import openmdao.api as om
+from mphys.multipoint import Multipoint
+from dafoam.mphys import DAFoamBuilder, OptFuncs
+from mphys.scenario_aerodynamic import ScenarioAerodynamic
+from pygeo.mphys import OM_DVGEOCOMP
+from pygeo import geo_utils
+from dafoam import PYDAFOAM
+
+
+# =============================================================================
+# Input Parameters
+# =============================================================================
+
+# NOTE: Input Parameters Must Also be Specified in "0/include/freestreamconditions"
+
+U0 = 233.603
+p0 = 108988.0
+T0 = 255.56
+nuTilda0 = 4.5e-5
+CL_target = 0.7
+twist0 = 2.52517169
+A0 = 0.01
+alpha0= 2.31
+k0= 8.5
+omega0= 5.6e4
+epsilon0= 4.3e4
+# rho is used for normalizing CD and CL
+rho0 = p0 / T0 / 287
+
+
+def calcUAndDir(UMag, alpha1):
+    dragDir = [float(np.cos(alpha1 * np.pi / 180)), float(np.sin(alpha1 * np.pi / 180)), 0.0]
+    liftDir = [float(-np.sin(alpha1 * np.pi / 180)), float(np.cos(alpha1 * np.pi / 180)), 0.0]
+    inletU = [float(UMag * np.cos(alpha1 * np.pi / 180)), float(UMag * np.sin(alpha1 * np.pi / 180)), 0.0]
+    return inletU, dragDir, liftDir
+
+
+inletU, dragDir, liftDir = calcUAndDir(U0, alpha0)
+
+
+# Input parameters for DAFoam
+daOptions = {
+    "designSurfaces": ["wing"],
+    "solverName": "DAHisaFoam",
+#    "primalInitCondition": {"U": [Ux, Uy, 0.0], "p": p0, "T": T0},
+    "primalMinResTol": 1.0e-7,
+    "primalBC": {
+        "U0": {"variable": "U", "patches": ["inout"], "value": inletU},
+        "p0": {"variable": "p", "patches": ["inout"], "value": [p0]},
+        "T0": {"variable": "T", "patches": ["inout"], "value": [T0]},
+        "nuTilda0": {"variable": "nuTilda", "patches": ["inout"], "value": [nuTilda0]},
+        "k0": {"variable": "k", "patches": ["inout"], "value": [k0]},
+        "omega0": {"variable": "omega", "patches": ["inout"], "value": [omega0]},
+        "epsilon0": {"variable": "epsilon", "patches": ["inout"], "value": [epsilon0]},
+        "useWallFunction": False,
+    },
+    "function": {
+        "CD": {
+            "type": "force",
+            "source": "patchToFace",
+            "patches": ["wing"],
+            "directionMode": "fixedDirection",
+            "direction": [1.0, 0.0, 0.0],
+            "scale": 1.0 / (0.5 * U0 * U0 * A0 * rho0),
+        },
+        "CL": {
+            "type": "force",
+            "source": "patchToFace",
+            "patches": ["wing"],
+            "directionMode": "fixedDirection",
+            "direction": [0.0, 1.0, 0.0],
+            "scale": 1.0 / (0.5 * U0 * U0 * A0 * rho0),
+        },
+    },
+    "checkMeshThreshold": {"maxNonOrth": 70.0, "maxSkewness": 6.0, "maxAspectRatio": 5000.0},
+
+}
+
+
+DASolver = PYDAFOAM(options=daOptions, comm=MPI.COMM_WORLD)
+DASolver()
+
+#
+#funcs = {}
+#evalFuncs = ["CD", "CL"]
+#DASolver.evalFunctions(funcs, evalFuncs)
diff --git a/RAE2822_Airfoil_HISA/system/controlDict b/RAE2822_Airfoil_HISA/system/controlDict
new file mode 100644
index 0000000..9ee733e
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/system/controlDict
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*---------------------------------*\ 
+| ========                 |                                                 | 
+| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           | 
+|  \    /   O peration     | Version:  v1812                                 | 
+|   \  /    A nd           | Web:      www.OpenFOAM.com                      | 
+|    \/     M anipulation  |                                                 | 
+\*--------------------------------------------------------------------------*/ 
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+application     hisa;
+startFrom       latestTime;
+startTime       0;
+stopAt          endTime;
+endTime         8000;
+deltaT          1;
+writeControl    timeStep;
+writeInterval   8000;
+purgeWrite      0;
+writeFormat     ascii;
+writePrecision  16;
+writeCompression on;
+timeFormat      general;
+timePrecision   8;
+runTimeModifiable false;
+HiSAPrint 0;
+
+DebugSwitches
+{
+    SolverPerformance 0;
+}
diff --git a/RAE2822_Airfoil_HISA/system/createPatchDict b/RAE2822_Airfoil_HISA/system/createPatchDict
new file mode 100644
index 0000000..5e2c407
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/system/createPatchDict
@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createPatchDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+pointSync false;
+
+// Patches to create.
+patches
+(
+    {
+        // Name of new patch
+        name symmetry1;
+
+        // Dictionary to construct new patch from
+        patchInfo
+        {
+            type symmetry;
+        }
+
+        // How to construct: either from 'patches' or 'set'
+        constructFrom patches;
+        patches (auto0);
+    }
+    {
+        // Name of new patch
+        name symmetry2;
+
+        // Dictionary to construct new patch from
+        patchInfo
+        {
+            type symmetry;
+        }
+
+        // How to construct: either from 'patches' or 'set'
+        constructFrom patches;
+        patches (auto1);
+    }
+
+    {
+        // Name of new patch
+        name wing;
+
+        // Dictionary to construct new patch from
+        patchInfo
+        {
+            type wall;
+        }
+
+        // How to construct: either from 'patches' or 'set'
+        constructFrom patches;
+        patches (auto2 auto3);
+    }
+
+    {
+        // Name of new patch
+        name inout;
+
+        // Dictionary to construct new patch from
+        patchInfo
+        {
+            type patch;
+        }
+
+        // How to construct: either from 'patches' or 'set'
+        constructFrom patches;
+        patches (auto4);
+    }
+);
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/system/decomposeParDict b/RAE2822_Airfoil_HISA/system/decomposeParDict
new file mode 100644
index 0000000..3f04af2
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/system/decomposeParDict
@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*---------------------------------*\ 
+| ========                 |                                                 | 
+| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           | 
+|  \    /   O peration     | Version:  v1812                                 | 
+|   \  /    A nd           | Web:      www.OpenFOAM.com                      | 
+|    \/     M anipulation  |                                                 | 
+\*--------------------------------------------------------------------------*/ 
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains     64;
+
+method                 scotch;
+
+simpleCoeffs 
+{ 
+    n                  (2 2 1);
+    delta              0.001;
+} 
+
+distributed            false;
+
+roots();
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/RAE2822_Airfoil_HISA/system/fvSchemes b/RAE2822_Airfoil_HISA/system/fvSchemes
new file mode 100644
index 0000000..f39c0a3
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/system/fvSchemes
@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+|                                                                             |
+|    HiSA: High Speed Aerodynamic solver                                      |
+|    Copyright (C) 2014-2017 Johan Heyns - CSIR, South Africa                 |
+|    Copyright (C) 2014-2017 Oliver Oxtoby - CSIR, South Africa               |
+|                                                                             |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+fluxScheme           JST; // AUSMPlusUp;
+jst_k2  0.5;
+jst_k4  0.02;
+sensorName pressure;
+lowMachAusm          false;
+
+ddtSchemes
+{
+    default          bounded dualTime rPseudoDeltaT steadyState;
+}
+
+gradSchemes
+{
+    default          faceLeastSquares linear;
+    grad(nuTilda)    Gauss linear;
+}
+
+divSchemes
+{
+    default          none;
+    div(tauMC)        Gauss linear;
+    div(phi,nuTilda) bounded Gauss upwind;
+    div(pc) bounded Gauss upwind;
+}
+
+laplacianSchemes
+{
+    default                     Gauss linear corrected;
+   laplacian(muEff,U)          Gauss linear compact;
+   laplacian(alphaEff,e)       Gauss linear compact;
+}
+
+interpolate         wVanLeer;
+interpolationSchemes
+{
+    default          linear;
+   reconstruct(rho) $interpolate;
+   reconstruct(U)   $interpolate;
+   reconstruct(T)   $interpolate;
+}
+
+snGradSchemes
+{
+    default          corrected;
+}
+
+wallDist
+{
+    method meshWave;
+}
+
+// ************************************************************************* //
diff --git a/RAE2822_Airfoil_HISA/system/fvSolution b/RAE2822_Airfoil_HISA/system/fvSolution
new file mode 100644
index 0000000..8b72658
--- /dev/null
+++ b/RAE2822_Airfoil_HISA/system/fvSolution
@@ -0,0 +1,74 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+|                                                                             |
+|    HiSA: High Speed Aerodynamic solver                                      |
+|    Copyright (C) 2014-2017 Johan Heyns - CSIR, South Africa                 |
+|    Copyright (C) 2014-2017 Oliver Oxtoby - CSIR, South Africa               |
+|                                                                             |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    nuTilda
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       0;
+        relTol          1e-1;
+        minIter         1;
+    }
+
+    yPsi
+    {
+        solver          GAMG;
+        smoother        GaussSeidel;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+        tolerance       0;
+        relTol          0;
+    }
+}
+
+relaxationFactors
+{
+    equations
+    {
+        nuTilda         0.5;
+        yWall           0.7;
+    }
+}
+
+flowSolver
+{
+    solver            GMRES;
+    GMRES
+    {
+        inviscidJacobian  LaxFriedrich;
+        viscousJacobian   laplacian;
+        preconditioner    LUSGS;
+
+        maxIter           20;
+        nKrylov           8;
+        solverTolRel      1e-1 (1e-1 1e-1 1e-1) 1e-1;
+    }
+}
+
+pseudoTime
+{
+    pseudoTol          1e-8 (1e-8 1e-8 1e-8) 1e-8;
+    pseudoCoNum        1.0;
+    pseudoCoNumMax     10.0;
+    localTimestepping  true;
+}
+
+// ************************************************************************* //