diff --git a/CHANGELOG.md b/CHANGELOG.md
index 40ba79164d2602656fcacb50cf50690060420ed2..39b31739797ecbfd2ccc783ccb22937e436b73d4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,10 @@ Changelog
 
 v4.4.0
 ------
+**Library:**
+- Add support for variable polynomial order for 3D simulations with continuous
+  Galerkin discretisation (!604)
+
 **IncNavierStokesSolver:**
 - Add ability to simulate additional scalar fields (!624)
 
diff --git a/library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp b/library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp
index 71cff2ba2128532595547f225388dad577a0e6cd..62e71867aab01003ba0f0b2928fe8a8406d1dc34 100644
--- a/library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp
+++ b/library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp
@@ -38,6 +38,7 @@
 #include <LocalRegions/Expansion.h>
 #include <LocalRegions/Expansion2D.h>
 #include <LocalRegions/Expansion3D.h>
+#include <LibUtilities/BasicUtils/ShapeType.hpp>
 
 
 #include <boost/config.hpp>
@@ -832,8 +833,8 @@ namespace Nektar
 
                 for(j = 0; j < nEdges; ++j)
                 {
-                    nEdgeIntCoeffs = exp->GetEdgeNcoeffs(j) - 2;
                     meshEdgeId = exp->GetGeom()->GetEid(j);
+                    nEdgeIntCoeffs = EdgeSize[meshEdgeId];
                     if(graph[1].count(meshEdgeId) == 0)
                     {
                         if(tempGraph[1].count(meshEdgeId) == 0)
@@ -1301,12 +1302,12 @@ namespace Nektar
             : AssemblyMap(pSession, variable)
         {
             int i, j, k, l;
+            int p, q, numModes0, numModes1;
             int cnt = 0;
             int intDofCnt;
             int meshVertId, meshEdgeId, meshFaceId;
             int globalId;
             int nEdgeInteriorCoeffs;
-            int nFaceInteriorCoeffs;
             int firstNonDirGraphVertId;
             LibUtilities::CommSharedPtr         vComm = m_comm->GetRowComm();
             LocalRegions::ExpansionSharedPtr    exp, bndExp;
@@ -1326,6 +1327,8 @@ namespace Nektar
             // Stores vertex, edge and face reordered vertices.
             DofGraph graph(3);
             DofGraph dofs(3);
+            vector<map<int, int> > faceModes(2);
+            map<int, LibUtilities::ShapeType> faceType;
 
             set<int> extraDirVerts, extraDirEdges;
             BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
@@ -1336,25 +1339,26 @@ namespace Nektar
             {
                 exp = locExpVector[i];
 
-                for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
+                for(j = 0; j < exp->GetNverts(); ++j)
                 {
                     dofs[0][exp->GetGeom()->GetVid(j)] = 1;
                 }
 
-                for(j = 0; j < locExpVector[i]->GetNedges(); ++j)
+                for(j = 0; j < exp->GetNedges(); ++j)
                 {
                     if (dofs[1].count(exp->GetGeom()->GetEid(j)) > 0)
                     {
                         if (dofs[1][exp->GetGeom()->GetEid(j)] !=
-                                    locExpVector[i]->GetEdgeNcoeffs(j)-2)
+                                    exp->GetEdgeNcoeffs(j)-2)
                         {
                             ASSERTL0( (exp->GetEdgeBasisType(j) == LibUtilities::eModified_A) ||
-                                      (exp->GetEdgeBasisType(j) == LibUtilities::eModified_B),
+                                      (exp->GetEdgeBasisType(j) == LibUtilities::eModified_B) ||
+                                      (exp->GetEdgeBasisType(j) == LibUtilities::eModified_C),
                                     "CG with variable order only available with modal expansion");
                         }
                         dofs[1][exp->GetGeom()->GetEid(j)] =
                                 min(dofs[1][exp->GetGeom()->GetEid(j)],
-                                    locExpVector[i]->GetEdgeNcoeffs(j)-2);
+                                    exp->GetEdgeNcoeffs(j)-2);
                     }
                     else
                     {
@@ -1363,30 +1367,37 @@ namespace Nektar
                     }
                 }
 
-                for(j = 0; j < locExpVector[i]->GetNfaces(); ++j)
+                for(j = 0; j < exp->GetNfaces(); ++j)
                 {
-                    if (dofs[2].count(exp->GetGeom()->GetFid(j)) > 0)
+                    faceOrient  = exp->GetGeom()->GetForient(j);
+                    meshFaceId  = exp->GetGeom()->GetFid(j);
+                    exp->GetFaceNumModes(j, faceOrient, numModes0, numModes1);
+
+                    if (faceModes[0].count(meshFaceId) > 0)
                     {
-                        if (dofs[2][exp->GetGeom()->GetFid(j)] !=
-                                     exp->GetFaceIntNcoeffs(j))
-                        {
-                            ASSERTL0( false,
-                                    "CG with variable order not available in 3D");
-                        }
-                        dofs[2][exp->GetGeom()->GetFid(j)] =
-                                min(dofs[2][exp->GetGeom()->GetFid(j)],
-                                    exp->GetFaceIntNcoeffs(j));
+                        faceModes[0][meshFaceId] =
+                                min(faceModes[0][meshFaceId], numModes0);
+
+                        faceModes[1][meshFaceId] =
+                                min(faceModes[1][meshFaceId], numModes1);
                     }
                     else
                     {
-                        dofs[2][exp->GetGeom()->GetFid(j)] =
-                            exp->GetFaceIntNcoeffs(j);
+                        faceModes[0][meshFaceId] = numModes0;
+                        faceModes[1][meshFaceId] = numModes1;
+
+                        // Get shape of this face
+                        SpatialDomains::Geometry3DSharedPtr geom;
+                        geom = boost::dynamic_pointer_cast<SpatialDomains::
+                                Geometry3D> (exp->GetGeom());
+                        faceType[meshFaceId] =
+                                geom->GetFace(j)->GetShapeType();
                     }
                 }
             }
             // Now use information from all partitions to determine
             //    the correct size
-            map<int, int>::iterator dofIt;
+            map<int, int>::iterator dofIt, dofIt2;
             // edges
             Array<OneD, long> edgeId (dofs[1].size());
             Array<OneD, NekDouble> edgeDof (dofs[1].size());
@@ -1403,19 +1414,41 @@ namespace Nektar
                 dofs[1][edgeId[i]] = (int) (edgeDof[i]+0.5);
             }
             // faces
-            Array<OneD, long> faceId (dofs[2].size());
-            Array<OneD, NekDouble> faceDof (dofs[2].size());
-            for(dofIt = dofs[2].begin(), i=0; dofIt != dofs[2].end(); dofIt++, i++)
+            Array<OneD, long> faceId (faceModes[0].size());
+            Array<OneD, NekDouble> faceP (faceModes[0].size());
+            Array<OneD, NekDouble> faceQ (faceModes[0].size());
+            for(dofIt = faceModes[0].begin(), dofIt2 = faceModes[1].begin(),i=0;
+                dofIt != faceModes[0].end(); dofIt++, dofIt2++, i++)
             {
                 faceId[i] = dofIt->first;
-                faceDof[i] = (NekDouble) dofIt->second;
+                faceP[i] = (NekDouble) dofIt->second;
+                faceQ[i] = (NekDouble) dofIt2->second;
             }
             Gs::gs_data *tmp2 = Gs::Init(faceId, vComm);
-            Gs::Gather(faceDof, Gs::gs_min, tmp2);
+            Gs::Gather(faceP, Gs::gs_min, tmp2);
+            Gs::Gather(faceQ, Gs::gs_min, tmp2);
             Gs::Finalise(tmp2);
-            for (i=0; i < dofs[2].size(); i++)
+            int P, Q;
+            for (i=0; i < faceModes[0].size(); i++)
             {
-                dofs[2][faceId[i]] = (int) (faceDof[i]+0.5);
+                faceModes[0][faceId[i]] = (int) (faceP[i]+0.5);
+                faceModes[1][faceId[i]] = (int) (faceQ[i]+0.5);
+                P = faceModes[0][faceId[i]];
+                Q = faceModes[1][faceId[i]];
+                if (faceType[faceId[i]] == LibUtilities::eQuadrilateral)
+                {
+                    // Quad face
+                    dofs[2][faceId[i]] =
+                      LibUtilities::StdQuadData::getNumberOfCoefficients(P,Q) -
+                      LibUtilities::StdQuadData::getNumberOfBndCoefficients(P,Q);
+                }
+                else
+                {
+                    // Tri face
+                    dofs[2][faceId[i]] =
+                      LibUtilities::StdTriData::getNumberOfCoefficients(P,Q) -
+                      LibUtilities::StdTriData::getNumberOfBndCoefficients(P,Q);
+                }
             }
 
             Array<OneD, const BndCond> bndCondVec(1, bndConditions);
@@ -1478,7 +1511,6 @@ namespace Nektar
 
                 for(j = 0; j < exp->GetNfaces(); ++j)
                 {
-                    nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
                     meshFaceId = exp->GetGeom()->GetFid(j);
                     graphVertOffset[graph[2][meshFaceId]+1] = dofs[2][meshFaceId];
                 }
@@ -1590,7 +1622,6 @@ namespace Nektar
 
                 for(j = 0; j < exp->GetNfaces(); ++j)
                 {
-                    nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
                     faceOrient          = exp->GetGeom()->GetForient(j);
                     meshFaceId          = exp->GetGeom()->GetFid(j);
 
@@ -1605,29 +1636,81 @@ namespace Nektar
                     exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
 
                     // Set the global DOF's for the interior modes of face j
-                    for(k = 0; k < dofs[2][exp->GetGeom()->GetFid(j)]; ++k)
-                    {
-                        m_localToGlobalMap[cnt+faceInteriorMap[k]] =
-                            graphVertOffset[graph[2][meshFaceId]]+k;
-                    }
-                    for(k = dofs[2][exp->GetGeom()->GetFid(j)]; k < nFaceInteriorCoeffs; ++k)
+                    exp->GetFaceNumModes(j, faceOrient, numModes0, numModes1);
+                    switch(faceType[meshFaceId])
                     {
-                        m_localToGlobalMap[cnt+faceInteriorMap[k]] =
-                            graphVertOffset[graph[2][meshFaceId]];
-                    }
-
-                    if(m_signChange)
+                    case LibUtilities::eQuadrilateral:
                     {
-                        for(k = 0; k < dofs[2][exp->GetGeom()->GetFid(j)]; ++k)
+                        int kLoc=0;
+                        k = 0;
+                        for( q = 2; q < numModes1; q++)
                         {
-                            m_localToGlobalSign[cnt+faceInteriorMap[k]] = (NekDouble) faceInteriorSign[k];
+                            for( p = 2; p < numModes0; p++)
+                            {
+                                if( (p < faceModes[0][meshFaceId]) &&
+                                    (q < faceModes[1][meshFaceId]))
+                                {
+                                    m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
+                                        graphVertOffset[graph[2][meshFaceId]]+k;
+                                    if(m_signChange)
+                                    {
+                                        m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] =
+                                            (NekDouble) faceInteriorSign[kLoc];
+                                    }
+                                    k++;
+                                }
+                                else
+                                {
+                                    m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
+                                        graphVertOffset[graph[2][meshFaceId]];
+                                    if(m_signChange)
+                                    {
+                                        m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] = 0.0;
+                                    }
+                                }
+                                kLoc++;
+                            }
                         }
-                        for(k = dofs[2][exp->GetGeom()->GetFid(j)]; k < nFaceInteriorCoeffs; ++k)
+                    }
+                    break;
+                    case LibUtilities::eTriangle:
+                    {
+                        int kLoc=0;
+                        k = 0;
+                        for( p = 2; p < numModes0; p++)
                         {
-                            m_localToGlobalSign[cnt+faceInteriorMap[k]] = 0;
+                            for( q = 1; q < numModes1-p; q++)
+                            {
+                                if( (p   < faceModes[0][meshFaceId]) &&
+                                    (p+q < faceModes[1][meshFaceId]))
+                                {
+                                    m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
+                                        graphVertOffset[graph[2][meshFaceId]]+k;
+                                    if(m_signChange)
+                                    {
+                                        m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] =
+                                            (NekDouble) faceInteriorSign[kLoc];
+                                    }
+                                    k++;
+                                }
+                                else
+                                {
+                                    m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
+                                        graphVertOffset[graph[2][meshFaceId]];
+                                    if(m_signChange)
+                                    {
+                                        m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] = 0.0;
+                                    }
+                                }
+                                kLoc++;
+                            }
                         }
                     }
-
+                    break;
+                    default:
+                        ASSERTL0(false,"Shape not recognised");
+                        break;
+                    }
                 }
             }
 
@@ -1688,19 +1771,28 @@ namespace Nektar
                         bndExp->GetEdgeInteriorMap(
                             k,edgeOrient,edgeInteriorMap,edgeInteriorSign);
 
-                        for(l = 0; l < nEdgeInteriorCoeffs; ++l)
+                        for(l = 0; l < dofs[1][meshEdgeId]; ++l)
                         {
                             m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
                                 graphVertOffset[graph[1][meshEdgeId]]+l;
                         }
+                        for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
+                        {
+                            m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
+                                graphVertOffset[graph[1][meshEdgeId]];
+                        }
 
                         // Fill the sign vector if required
                         if(m_signChange)
                         {
-                            for(l = 0; l < nEdgeInteriorCoeffs; ++l)
+                            for(l = 0; l < dofs[1][meshEdgeId]; ++l)
                             {
                                 m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = (NekDouble) edgeInteriorSign[l];
                             }
+                            for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
+                            {
+                                m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = 0.0;
+                            }
                         }
 
                         if (bndConditions[i]->GetBoundaryConditionType() !=
@@ -1714,7 +1806,7 @@ namespace Nektar
                             foundExtraEdges.count(meshEdgeId) == 0 &&
                             nEdgeInteriorCoeffs > 0)
                         {
-                            for(l = 0; l < nEdgeInteriorCoeffs; ++l)
+                            for(l = 0; l < dofs[1][meshEdgeId]; ++l)
                             {
                                 int loc = bndCondExp[i]->GetCoeff_Offset(j) +
                                     edgeInteriorMap[l];
@@ -1723,6 +1815,15 @@ namespace Nektar
                                 ExtraDirDof t(loc, gid, edgeInteriorSign[l]);
                                 m_extraDirDofs[i].push_back(t);
                             }
+                            for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
+                            {
+                                int loc = bndCondExp[i]->GetCoeff_Offset(j) +
+                                    edgeInteriorMap[l];
+                                int gid = graphVertOffset[
+                                    graph[1][meshEdgeId]];
+                                ExtraDirDof t(loc, gid, 0.0);
+                                m_extraDirDofs[i].push_back(t);
+                            }
                             foundExtraEdges.insert(meshEdgeId);
                         }
                     }
@@ -2007,7 +2108,7 @@ namespace Nektar
             int maxIntDof = 0;
             int dof = 0;
             int cnt;
-            int i,j,k;
+            int i,j,k,l;
             int meshVertId;
             int meshEdgeId;
             int meshFaceId;
@@ -2108,14 +2209,22 @@ namespace Nektar
                     dof = exp->GetEdgeNcoeffs(j)-2;
 
                     // Set the global DOF's for the interior modes of edge j
-                    //    run backwards because of variable P case "ghost" modes
-                    for(k = dof-1; k >= 0; --k)
+                    //    for varP, ignore modes with sign == 0
+                    for(k = 0, l = 0; k < dof; ++k)
                     {
+                        if (m_signChange)
+                        {
+                            if (m_localToGlobalSign[cnt+edgeInteriorMap[k]]==0)
+                            {
+                                continue;
+                            }
+                        }
                         vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
                         m_globalToUniversalMap[vGlobalId]
-                           = nVert + meshEdgeId * maxEdgeDof + k + 1;
+                           = nVert + meshEdgeId * maxEdgeDof + l + 1;
                         m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
                         maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
+                        l++;
                     }
                 }
 
@@ -2140,16 +2249,23 @@ namespace Nektar
                     exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
                     dof = exp->GetFaceIntNcoeffs(j);
 
-
-                    for(k = dof-1; k >= 0; --k)
+                    for(k = 0, l = 0; k < dof; ++k)
                     {
+                        if (m_signChange)
+                        {
+                            if (m_localToGlobalSign[cnt+faceInteriorMap[k]]==0)
+                            {
+                                continue;
+                            }
+                        }
                         vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
                         m_globalToUniversalMap[vGlobalId]
                            = nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
-                                   + k + 1;
+                                   + l + 1;
                         m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
 
                         maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
+                        l++;
                     }
                 }
 
diff --git a/library/MultiRegions/ContField3D.cpp b/library/MultiRegions/ContField3D.cpp
index 556f775688bd20d0fee8c0404d401b5738d6f05b..9dfd300af9b946ce5b986280ad3949af5780670b 100644
--- a/library/MultiRegions/ContField3D.cpp
+++ b/library/MultiRegions/ContField3D.cpp
@@ -339,7 +339,14 @@ namespace Nektar
                   for(int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
                   {
                       sign = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
-                      inout[map[bndcnt++]] = sign * coeffs[j];
+                      if(sign)
+                      {
+                          inout[map[bndcnt++]] = sign * coeffs[j];
+                      }
+                      else
+                      {
+                          bndcnt++;
+                      }
                   }
               }
               else
@@ -465,7 +472,14 @@ namespace Nektar
                   {
                       sign = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsSign(
                           bndcnt);
-                      tmp[bndMap[bndcnt++]] = sign * coeffs[j];
+                      if (sign)
+                      {
+                          tmp[bndMap[bndcnt++]] = sign * coeffs[j];
+                      }
+                      else
+                      {
+                          bndcnt++;
+                      }
                   }
               }
               else
diff --git a/library/SolverUtils/DriverAdaptive.cpp b/library/SolverUtils/DriverAdaptive.cpp
index 8b922ef16eb50e4492c202ba4d4c521910b82e99..b7d42e3dd7d04b93741004c56c020b1ea3ec4c57 100644
--- a/library/SolverUtils/DriverAdaptive.cpp
+++ b/library/SolverUtils/DriverAdaptive.cpp
@@ -38,6 +38,10 @@
 #include <SolverUtils/DriverAdaptive.h>
 #include <StdRegions/StdQuadExp.h>
 #include <StdRegions/StdTriExp.h>
+#include <StdRegions/StdTetExp.h>
+#include <StdRegions/StdPrismExp.h>
+#include <StdRegions/StdPyrExp.h>
+#include <StdRegions/StdHexExp.h>
 #include <GlobalMapping/Mapping.h>
 
 using namespace std;
@@ -115,6 +119,7 @@ void DriverAdaptive::v_Execute(ostream &out)
         nExp    = m_equ[0]->UpdateFields()[0]->GetExpSize();
         nPlanes = 1;
     }
+    int  expdim   = m_equ[0]->UpdateFields()[0]->GetGraph()->GetMeshDimension();
 
     int       nFields  = m_equ[0]->UpdateFields().num_elements();
     int       numSteps = m_session->GetParameter("NumSteps");
@@ -128,6 +133,9 @@ void DriverAdaptive::v_Execute(ostream &out)
                                            m_equ[0]->UpdateFields());
 
     // Adaptive loop
+    Array<OneD, int> P(expdim);
+    Array<OneD, int> numPoints(expdim);
+    Array<OneD, LibUtilities::PointsKey> ptsKey(expdim);
     LocalRegions::ExpansionSharedPtr Exp;
     for (int i = 1; i < nRuns; i++)
     {
@@ -143,39 +151,91 @@ void DriverAdaptive::v_Execute(ostream &out)
             offset = fields[sensorVar]->GetPhys_Offset(n);
             Exp    = fields[sensorVar]->GetExp(n);
 
-            int P  = Exp->GetBasis(0)->GetNumModes();
-            int Q  = Exp->GetBasis(1)->GetNumModes();
-            int qa = Exp->GetBasis(0)->GetNumPoints();
-            int qb = Exp->GetBasis(1)->GetNumPoints();
+            for( int k = 0; k < expdim; ++k)
+            {
+                P[k]         = Exp->GetBasis(k)->GetNumModes();
+                numPoints[k] = Exp->GetBasis(k)->GetNumPoints();
+                ptsKey[k]    = LibUtilities::PointsKey (numPoints[k],
+                               Exp->GetBasis(k)->GetPointsType());
+            }
 
             // Declare orthogonal basis.
-            LibUtilities::PointsKey pa(qa, Exp->GetBasis(0)->GetPointsType());
-            LibUtilities::PointsKey pb(qb, Exp->GetBasis(1)->GetPointsType());
-
             StdRegions::StdExpansionSharedPtr OrthoExp;
             switch (Exp->GetGeom()->GetShapeType())
             {
                 case LibUtilities::eQuadrilateral:
                 {
-                    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, P - 1,
-                                              pa);
-                    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, Q - 1,
-                                              pb);
+                    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, P[0] - 1,
+                                              ptsKey[0]);
+                    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, P[1] - 1,
+                                              ptsKey[1]);
                     OrthoExp = MemoryManager<
                         StdRegions::StdQuadExp>::AllocateSharedPtr(Ba, Bb);
                     break;
                 }
                 case LibUtilities::eTriangle:
                 {
-                    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, P - 1,
-                                              pa);
-                    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B, Q - 1,
-                                              pb);
+                    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, P[0] - 1,
+                                              ptsKey[0]);
+                    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B, P[1] - 1,
+                                              ptsKey[1]);
                     OrthoExp =
                         MemoryManager<StdRegions::StdTriExp>::AllocateSharedPtr(
                             Ba, Bb);
                     break;
                 }
+                case LibUtilities::eTetrahedron:
+                {
+                    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, P[0] - 1,
+                                              ptsKey[0]);
+                    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B, P[1] - 1,
+                                              ptsKey[1]);
+                    LibUtilities::BasisKey Bc(LibUtilities::eOrtho_C, P[2] - 1,
+                                              ptsKey[2]);
+                    OrthoExp =
+                        MemoryManager<StdRegions::StdTetExp>::AllocateSharedPtr(
+                            Ba, Bb, Bc);
+                    break;
+                }
+                case LibUtilities::ePyramid:
+                {
+                    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, P[0] - 1,
+                                              ptsKey[0]);
+                    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, P[1] - 1,
+                                              ptsKey[1]);
+                    LibUtilities::BasisKey Bc(LibUtilities::eOrtho_C, P[2] - 1,
+                                              ptsKey[2]);
+                    OrthoExp =
+                        MemoryManager<StdRegions::StdPyrExp>::AllocateSharedPtr(
+                            Ba, Bb, Bc);
+                    break;
+                }
+                case LibUtilities::ePrism:
+                {
+                    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, P[0] - 1,
+                                              ptsKey[0]);
+                    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, P[1] - 1,
+                                              ptsKey[1]);
+                    LibUtilities::BasisKey Bc(LibUtilities::eOrtho_B, P[2] - 1,
+                                              ptsKey[2]);
+                    OrthoExp =
+                        MemoryManager<StdRegions::StdPrismExp>::AllocateSharedPtr(
+                            Ba, Bb, Bc);
+                    break;
+                }
+                case LibUtilities::eHexahedron:
+                {
+                    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, P[0] - 1,
+                                              ptsKey[0]);
+                    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, P[1] - 1,
+                                              ptsKey[1]);
+                    LibUtilities::BasisKey Bc(LibUtilities::eOrtho_A, P[2] - 1,
+                                              ptsKey[2]);
+                    OrthoExp =
+                        MemoryManager<StdRegions::StdHexExp>::AllocateSharedPtr(
+                            Ba, Bb, Bc);
+                    break;
+                }
                 default:
                     ASSERTL0(false, "Shape not supported.");
                     break;
@@ -272,11 +332,11 @@ void DriverAdaptive::v_Execute(ostream &out)
             }
 
             // Determine what to do with the polynomial order
-            if ((error > upperTol) && (P < maxP))
+            if ((error > upperTol) && (P[0] < maxP))
             {
                 deltaP[Exp->GetGeom()->GetGlobalID()] = 1;
             }
-            else if ((error < lowerTol) && P > minP)
+            else if ((error < lowerTol) && P[0] > minP)
             {
                 deltaP[Exp->GetGeom()->GetGlobalID()] = -1;
             }
@@ -387,20 +447,19 @@ void DriverAdaptive::ReplaceExpansion(
     map<int, int>                                deltaP)
 {
     int nExp, nDim;
+    int  expdim   = m_equ[0]->UpdateFields()[0]->GetGraph()->GetMeshDimension();
 
     // Get field definitions
     std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddefs =
         fields[0]->GetFieldDefinitions();
 
-    int expDim = 2;
-
     if (fielddefs[0]->m_numHomogeneousDir == 1)
     {
-        nDim = 3;
+        nDim = expdim+1;
     }
     else
     {
-        nDim = 2;
+        nDim = expdim;
     }
 
     // Add variables to field definition
@@ -569,7 +628,7 @@ void DriverAdaptive::ReplaceExpansion(
             {
                 eId = fielddefs[f]->m_elementIDs[n];
 
-                for (int i = 0; i < expDim; i++)
+                for (int i = 0; i < expdim; i++)
                 {
                     order[i] = deltaP[eId];
                 }
diff --git a/library/StdRegions/StdExpansion.cpp b/library/StdRegions/StdExpansion.cpp
index 322946a5b7554d7a5e8eb5621cda862e6e53f130..2bcd02c9774fdd4081ed84cd656c9078f9d734b9 100644
--- a/library/StdRegions/StdExpansion.cpp
+++ b/library/StdRegions/StdExpansion.cpp
@@ -1414,7 +1414,7 @@ namespace Nektar
                 NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
             }
 
-        int StdExpansion::v_GetVertexMap(const int localVertexId,
+            int StdExpansion::v_GetVertexMap(const int localVertexId,
                                          bool useCoeffPacking)
             {
                 NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
@@ -1428,6 +1428,15 @@ namespace Nektar
                 NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
             }
 
+            void StdExpansion::v_GetFaceNumModes(
+                                              const int fid,
+                                              const Orientation faceOrient,
+                                              int &numModes0,
+                                              int &numModes1)
+            {
+                NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
+            }
+
             void StdExpansion::v_GetFaceInteriorMap(const int fid, const Orientation faceOrient,
                                               Array<OneD, unsigned int> &maparray,
                                               Array<OneD, int> &signarray)
diff --git a/library/StdRegions/StdExpansion.h b/library/StdRegions/StdExpansion.h
index 9e12fc32a84a9d91564dd81ba1287c27e4ff0cc3..f23a6bb46d62a2e85b7fd6f87d991d5492849ab9 100644
--- a/library/StdRegions/StdExpansion.h
+++ b/library/StdRegions/StdExpansion.h
@@ -836,6 +836,13 @@ namespace Nektar
                 v_GetEdgeInteriorMap(eid,edgeOrient,maparray,signarray);
             }
 
+            void GetFaceNumModes(const int fid, const Orientation faceOrient,
+                                    int &numModes0,
+                                    int &numModes1)
+            {
+                v_GetFaceNumModes(fid,faceOrient,numModes0,numModes1);
+            }
+
             void GetFaceInteriorMap(const int fid, const Orientation faceOrient,
                                     Array<OneD, unsigned int> &maparray,
                                     Array<OneD, int> &signarray)
@@ -1706,6 +1713,12 @@ namespace Nektar
                                               Array<OneD, unsigned int> &maparray,
                                               Array<OneD, int> &signarray);
 
+            STD_REGIONS_EXPORT virtual void v_GetFaceNumModes(
+                                              const int fid,
+                                              const Orientation faceOrient,
+                                              int &numModes0,
+                                              int &numModes1);
+
             STD_REGIONS_EXPORT virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient,
                                               Array<OneD, unsigned int> &maparray,
                                               Array<OneD, int> &signarray);
diff --git a/library/StdRegions/StdHexExp.cpp b/library/StdRegions/StdHexExp.cpp
index 92ad7df22b9907cfa8ba149bc11e57209c527010..44338018c1a85dd75d154bf99f221e898fa9f663 100644
--- a/library/StdRegions/StdHexExp.cpp
+++ b/library/StdRegions/StdHexExp.cpp
@@ -927,6 +927,50 @@ namespace Nektar
             }
         }
 
+        void StdHexExp::v_GetFaceNumModes(
+                    const int                  fid,
+                    const Orientation          faceOrient,
+                    int &numModes0,
+                    int &numModes1)
+        {
+            int nummodes [3] = {m_base[0]->GetNumModes(),
+                                m_base[1]->GetNumModes(),
+                                m_base[2]->GetNumModes()};
+            switch(fid)
+            {
+            case 0:
+            case 5:
+                {
+                    numModes0 = nummodes[0];
+                    numModes1 = nummodes[1];
+                }
+                break;
+            case 1:
+            case 3:
+                {
+                    numModes0 = nummodes[0];
+                    numModes1 = nummodes[2];
+                }
+                break;
+            case 2:
+            case 4:
+                {
+                    numModes0 = nummodes[1];
+                    numModes1 = nummodes[2];
+                }
+                break;
+            default:
+                {
+                    ASSERTL0(false,"fid out of range");
+                }
+                break;
+            }
+
+            if ( faceOrient >= 9 )
+            {
+                std::swap(numModes0, numModes1);
+            }
+        }
 
         /**
          * Only for basis type Modified_A or GLL_LAGRANGE in all directions.
diff --git a/library/StdRegions/StdHexExp.h b/library/StdRegions/StdHexExp.h
index df3baf52a65f082d4eab6fdbfaa062b378ca7f1d..8e31434548a758a855a76bc37f8e0c92af4bfdd9 100644
--- a/library/StdRegions/StdHexExp.h
+++ b/library/StdRegions/StdHexExp.h
@@ -199,7 +199,11 @@ namespace Nektar
                     Array<OneD, NekDouble> &coords_x,
                     Array<OneD, NekDouble> &coords_y,
                     Array<OneD, NekDouble> &coords_z);
-
+            STD_REGIONS_EXPORT virtual void v_GetFaceNumModes(
+                    const int                  fid,
+                    const Orientation          faceOrient,
+                    int &numModes0,
+                    int &numModes1);
 
             //--------------------------
             // Mappings
diff --git a/library/StdRegions/StdPrismExp.cpp b/library/StdRegions/StdPrismExp.cpp
index 90f19e3f0e033d34eceba3555c21908e37b31c97..a05328731a93c43342efc8cfd3c8e0fff70dfdde 100644
--- a/library/StdRegions/StdPrismExp.cpp
+++ b/library/StdRegions/StdPrismExp.cpp
@@ -728,6 +728,47 @@ namespace Nektar
             StdPrismExp::v_BwdTrans(tmp, outarray);
         }
 
+        void StdPrismExp::v_GetFaceNumModes(
+                    const int                  fid,
+                    const Orientation          faceOrient,
+                    int &numModes0,
+                    int &numModes1)
+        {
+            int nummodes [3] = {m_base[0]->GetNumModes(),
+                                m_base[1]->GetNumModes(),
+                                m_base[2]->GetNumModes()};
+            switch(fid)
+            {
+            // base quad
+            case 0:
+                {
+                    numModes0 = nummodes[0];
+                    numModes1 = nummodes[1];
+                }
+                break;
+            // front and back quad
+            case 2:
+            case 4:
+                {
+                    numModes0 = nummodes[1];
+                    numModes1 = nummodes[2];
+                }
+                break;
+            // triangles
+            case 1:
+            case 3:
+                {
+                    numModes0 = nummodes[0];
+                    numModes1 = nummodes[2];
+                }
+                break;
+            }
+
+            if ( faceOrient >= 9 )
+            {
+                std::swap(numModes0, numModes1);
+            }
+        }
 
         //---------------------------------------
         // Helper functions
diff --git a/library/StdRegions/StdPrismExp.h b/library/StdRegions/StdPrismExp.h
index 9b82bc8993a40bf4065d4625b1f67ee1e1b41f37..31ac7b45b78e084769d69d212fe5da8131a3ff38 100644
--- a/library/StdRegions/StdPrismExp.h
+++ b/library/StdRegions/StdPrismExp.h
@@ -165,6 +165,11 @@ namespace Nektar
             STD_REGIONS_EXPORT virtual void v_FillMode(
                 const int mode,
                 Array<OneD, NekDouble> &outarray);
+            STD_REGIONS_EXPORT virtual void v_GetFaceNumModes(
+                    const int                  fid,
+                    const Orientation          faceOrient,
+                    int &numModes0,
+                    int &numModes1);
 
             //---------------------------------------
             // Helper functions
diff --git a/library/StdRegions/StdPyrExp.cpp b/library/StdRegions/StdPyrExp.cpp
index b2325c51aeb9146cf0610ef66c51e551667a8e52..911eaed77a761156a7e9c3c21e64bd454ea50363 100644
--- a/library/StdRegions/StdPyrExp.cpp
+++ b/library/StdRegions/StdPyrExp.cpp
@@ -867,8 +867,47 @@ namespace Nektar
             tmp[mode] = 1.0;
             v_BwdTrans(tmp, outarray);
         }
-        
-        
+
+        void StdPyrExp::v_GetFaceNumModes(
+                    const int                  fid,
+                    const Orientation          faceOrient,
+                    int &numModes0,
+                    int &numModes1)
+        {
+            int nummodes [3] = {m_base[0]->GetNumModes(),
+                                m_base[1]->GetNumModes(),
+                                m_base[2]->GetNumModes()};
+            switch(fid)
+            {
+            // quad
+            case 0:
+                {
+                    numModes0 = nummodes[0];
+                    numModes1 = nummodes[1];
+                }
+                break;
+            case 1:
+            case 3:
+                {
+                    numModes0 = nummodes[0];
+                    numModes1 = nummodes[2];
+                }
+                break;
+            case 2:
+            case 4:
+                {
+                    numModes0 = nummodes[1];
+                    numModes1 = nummodes[2];
+                }
+                break;
+            }
+
+            if ( faceOrient >= 9 )
+            {
+                std::swap(numModes0, numModes1);
+            }
+        }
+
         //---------------------------------------
         // Helper functions
         //---------------------------------------
diff --git a/library/StdRegions/StdPyrExp.h b/library/StdRegions/StdPyrExp.h
index 63683694aa36641c6e91066927ca91ecc9b265f5..4b93af8efac1e659b1c9499b34448fbd79b83021 100644
--- a/library/StdRegions/StdPyrExp.h
+++ b/library/StdRegions/StdPyrExp.h
@@ -191,6 +191,11 @@ namespace Nektar
             STD_REGIONS_EXPORT virtual void v_FillMode(
                 const int mode, 
                 Array<OneD, NekDouble> &outarray);  
+            STD_REGIONS_EXPORT virtual void v_GetFaceNumModes(
+                    const int                  fid,
+                    const Orientation          faceOrient,
+                    int &numModes0,
+                    int &numModes1);
 
             //---------------------------------------
             // Helper functions
diff --git a/library/StdRegions/StdTetExp.cpp b/library/StdRegions/StdTetExp.cpp
index 8a402c7a3cd9b2244370e403f6e92c44f0e2ca32..98d5dce106529a78c81ee89d5bc8a8e7044d2296 100644
--- a/library/StdRegions/StdTetExp.cpp
+++ b/library/StdRegions/StdTetExp.cpp
@@ -906,6 +906,38 @@ namespace Nektar
             StdTetExp::v_BwdTrans(tmp, outarray);
         }
 
+        void StdTetExp::v_GetFaceNumModes(
+                    const int                  fid,
+                    const Orientation          faceOrient,
+                    int &numModes0,
+                    int &numModes1)
+        {
+            int nummodes [3] = {m_base[0]->GetNumModes(),
+                                m_base[1]->GetNumModes(),
+                                m_base[2]->GetNumModes()};
+            switch(fid)
+            {
+            case 0:
+                {
+                    numModes0 = nummodes[0];
+                    numModes1 = nummodes[1];
+                }
+                break;
+            case 1:
+                {
+                    numModes0 = nummodes[0];
+                    numModes1 = nummodes[2];
+                }
+                break;
+            case 2:
+            case 3:
+                {
+                    numModes0 = nummodes[1];
+                    numModes1 = nummodes[2];
+                }
+                break;
+            }
+        }
 
         //---------------------------
         // Helper functions
diff --git a/library/StdRegions/StdTetExp.h b/library/StdRegions/StdTetExp.h
index 224329c2e8ea9c85a362ae00e918d7204a6f236f..9c6d2d44955ff5070b181ad9c90b001be1e60206 100644
--- a/library/StdRegions/StdTetExp.h
+++ b/library/StdRegions/StdTetExp.h
@@ -173,7 +173,11 @@ namespace Nektar
             STD_REGIONS_EXPORT virtual void v_FillMode(
                 const int                     mode,
                       Array<OneD, NekDouble>& outarray);
-
+            STD_REGIONS_EXPORT virtual void v_GetFaceNumModes(
+                    const int                  fid,
+                    const Orientation          faceOrient,
+                    int &numModes0,
+                    int &numModes1);
 
             //---------------------------
             // Helper functions
diff --git a/solvers/IncNavierStokesSolver/CMakeLists.txt b/solvers/IncNavierStokesSolver/CMakeLists.txt
index 657474305f550f94e049c5a83666d0d64b8438f6..423da621a8895ba23919c3f164776c2f5045a05f 100644
--- a/solvers/IncNavierStokesSolver/CMakeLists.txt
+++ b/solvers/IncNavierStokesSolver/CMakeLists.txt
@@ -49,7 +49,9 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
     ADD_NEKTAR_TEST(Channel_Flow_3modes_rad)    
     ADD_NEKTAR_TEST(Couette_3DH2D_MVM)
     ADD_NEKTAR_TEST(Hex_channel_m3)
+    ADD_NEKTAR_TEST(Hex_channel_varP)
     ADD_NEKTAR_TEST(Pyr_channel_m3)
+    ADD_NEKTAR_TEST(Pyr_channel_varP)
     ADD_NEKTAR_TEST_LENGTHY(Hex_channel_m8)
     ADD_NEKTAR_TEST_LENGTHY(Hex_channel_m8_srhs)
     ADD_NEKTAR_TEST(Hex_Kovasnay_SubStep)
@@ -79,10 +81,12 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
     ADD_NEKTAR_TEST(Tet_channel_m3_dealiasing)
     ADD_NEKTAR_TEST(Tet_channel_robin_m3)
     ADD_NEKTAR_TEST(Tet_channel_m4_per)
+    ADD_NEKTAR_TEST(Tet_channel_varP)
     ADD_NEKTAR_TEST_LENGTHY(Tet_prism_channel_rad)
     ADD_NEKTAR_TEST_LENGTHY(Tet_equitri)
     ADD_NEKTAR_TEST_LENGTHY(Prism_channel_m6)
     ADD_NEKTAR_TEST_LENGTHY(Prism_channel_m6_dealiasing)
+    ADD_NEKTAR_TEST(Prism_channel_varP)
     ADD_NEKTAR_TEST(ChannelSpongeNSE)
     ADD_NEKTAR_TEST(ChannelSpongeLNSE)
     ADD_NEKTAR_TEST(ChanFlow_Standard_BodyForce)
@@ -119,6 +123,7 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
         ADD_NEKTAR_TEST(ChanFlow_m3_par)
         ADD_NEKTAR_TEST_LENGTHY(ChanFlow_m8_BodyForce_par)
         ADD_NEKTAR_TEST_LENGTHY(Hex_channel_m8_par)
+        ADD_NEKTAR_TEST(Hex_channel_varP_par)
         ADD_NEKTAR_TEST_LENGTHY(Pyr_channel_m6_par)
         ADD_NEKTAR_TEST_LENGTHY(Hex_channel_m8_srhs_par)
         ADD_NEKTAR_TEST(Tet_channel_m4_per_xxt_ml_par)
diff --git a/solvers/IncNavierStokesSolver/Tests/Hex_channel_varP.tst b/solvers/IncNavierStokesSolver/Tests/Hex_channel_varP.tst
new file mode 100644
index 0000000000000000000000000000000000000000..16f0efd017ca63752824cc07d3565f1e7f9c1e58
--- /dev/null
+++ b/solvers/IncNavierStokesSolver/Tests/Hex_channel_varP.tst
@@ -0,0 +1,24 @@
+<?xml version="1.0" encoding="utf-8"?>
+<test>
+    <description>3D channel flow, Hexahedral elements, variable P</description>
+    <executable>IncNavierStokesSolver</executable>
+    <parameters>Hex_channel_varP.xml</parameters>
+    <files>
+        <file description="Session File">Hex_channel_varP.xml</file>
+    </files>
+    <metrics>
+        <metric type="L2" id="1">
+            <value variable="u" tolerance="1e-12">0</value>
+            <value variable="v" tolerance="1e-12">0</value>
+            <value variable="w" tolerance="1e-12">0</value>
+            <value variable="p" tolerance="1e-11">0</value>
+        </metric>
+        <metric type="Linf" id="2">
+            <value variable="u" tolerance="1e-12">0</value>
+            <value variable="v" tolerance="1e-12">0</value>
+            <value variable="w" tolerance="1e-12">0</value>
+            <value variable="p" tolerance="1e-11">0</value>
+        </metric>
+    </metrics>
+</test>
+
diff --git a/solvers/IncNavierStokesSolver/Tests/Hex_channel_varP.xml b/solvers/IncNavierStokesSolver/Tests/Hex_channel_varP.xml
new file mode 100644
index 0000000000000000000000000000000000000000..04361fab01ab61581dc85dae14328f8c5f3e1864
--- /dev/null
+++ b/solvers/IncNavierStokesSolver/Tests/Hex_channel_varP.xml
@@ -0,0 +1,242 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <GEOMETRY DIM="3" SPACE="3">
+        <VERTEX>
+            <V ID="0">0.00000000e+00 0.00000000e+00 0.00000000e+00</V>
+            <V ID="1">5.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="2">5.00000000e-01 5.00000000e-01 0.00000000e+00</V>
+            <V ID="3">0.00000000e+00 5.00000000e-01 0.00000000e+00</V>
+            <V ID="4">0.00000000e+00 0.00000000e+00 5.00000000e-01</V>
+            <V ID="5">5.00000000e-01 0.00000000e+00 5.00000000e-01</V>
+            <V ID="6">5.00000000e-01 5.00000000e-01 5.00000000e-01</V>
+            <V ID="7">0.00000000e+00 5.00000000e-01 5.00000000e-01</V>
+            <V ID="8">0.00000000e+00 0.00000000e+00 1.00000000e+00</V>
+            <V ID="9">5.00000000e-01 0.00000000e+00 1.00000000e+00</V>
+            <V ID="10">5.00000000e-01 5.00000000e-01 1.00000000e+00</V>
+            <V ID="11">0.00000000e+00 5.00000000e-01 1.00000000e+00</V>
+            <V ID="12">5.00000000e-01 1.00000000e+00 0.00000000e+00</V>
+            <V ID="13">0.00000000e+00 1.00000000e+00 0.00000000e+00</V>
+            <V ID="14">5.00000000e-01 1.00000000e+00 5.00000000e-01</V>
+            <V ID="15">0.00000000e+00 1.00000000e+00 5.00000000e-01</V>
+            <V ID="16">5.00000000e-01 1.00000000e+00 1.00000000e+00</V>
+            <V ID="17">0.00000000e+00 1.00000000e+00 1.00000000e+00</V>
+            <V ID="18">1.00000000e+00 0.00000000e+00 0.00000000e+00</V>
+            <V ID="19">1.00000000e+00 5.00000000e-01 0.00000000e+00</V>
+            <V ID="20">1.00000000e+00 0.00000000e+00 5.00000000e-01</V>
+            <V ID="21">1.00000000e+00 5.00000000e-01 5.00000000e-01</V>
+            <V ID="22">1.00000000e+00 0.00000000e+00 1.00000000e+00</V>
+            <V ID="23">1.00000000e+00 5.00000000e-01 1.00000000e+00</V>
+            <V ID="24">1.00000000e+00 1.00000000e+00 0.00000000e+00</V>
+            <V ID="25">1.00000000e+00 1.00000000e+00 5.00000000e-01</V>
+            <V ID="26">1.00000000e+00 1.00000000e+00 1.00000000e+00</V>
+        </VERTEX>
+        <EDGE>
+            <E ID="0">    0  1   </E>
+            <E ID="1">    0  3   </E>
+            <E ID="2">    0  4   </E>
+            <E ID="3">    1  2   </E>
+            <E ID="4">    1  5   </E>
+            <E ID="5">    2  3   </E>
+            <E ID="6">    2  6   </E>
+            <E ID="7">    3  7   </E>
+            <E ID="8">    4  5   </E>
+            <E ID="9">    4  7   </E>
+            <E ID="10">    5  6   </E>
+            <E ID="11">    6  7   </E>
+            <E ID="12">    4  8   </E>
+            <E ID="13">    5  9   </E>
+            <E ID="14">    6  10   </E>
+            <E ID="15">    7  11   </E>
+            <E ID="16">    8  9   </E>
+            <E ID="17">    8  11   </E>
+            <E ID="18">    9  10   </E>
+            <E ID="19">   10  11   </E>
+            <E ID="20">    3  13   </E>
+            <E ID="21">    2  12   </E>
+            <E ID="22">   12  13   </E>
+            <E ID="23">   12  14   </E>
+            <E ID="24">   13  15   </E>
+            <E ID="25">    7  15   </E>
+            <E ID="26">    6  14   </E>
+            <E ID="27">   14  15   </E>
+            <E ID="28">   14  16   </E>
+            <E ID="29">   15  17   </E>
+            <E ID="30">   11  17   </E>
+            <E ID="31">   10  16   </E>
+            <E ID="32">   16  17   </E>
+            <E ID="33">    1  18   </E>
+            <E ID="34">   18  19   </E>
+            <E ID="35">   18  20   </E>
+            <E ID="36">   19  2   </E>
+            <E ID="37">   19  21   </E>
+            <E ID="38">    5  20   </E>
+            <E ID="39">   20  21   </E>
+            <E ID="40">   21  6   </E>
+            <E ID="41">   20  22   </E>
+            <E ID="42">   21  23   </E>
+            <E ID="43">    9  22   </E>
+            <E ID="44">   22  23   </E>
+            <E ID="45">   23  10   </E>
+            <E ID="46">   19  24   </E>
+            <E ID="47">   24  12   </E>
+            <E ID="48">   24  25   </E>
+            <E ID="49">   21  25   </E>
+            <E ID="50">   25  14   </E>
+            <E ID="51">   25  26   </E>
+            <E ID="52">   23  26   </E>
+            <E ID="53">   26  16   </E>
+        </EDGE>
+        <FACE>
+            <Q ID="0">         0         3         5         1</Q>
+            <Q ID="1">         0         4         8         2</Q>
+            <Q ID="2">         3         6        10         4</Q>
+            <Q ID="3">         5         7        11         6</Q>
+            <Q ID="4">         1         2         9         7</Q>
+            <Q ID="5">         8        10        11         9</Q>
+            <Q ID="6">         8        13        16        12</Q>
+            <Q ID="7">        10        14        18        13</Q>
+            <Q ID="8">        11        15        19        14</Q>
+            <Q ID="9">         9        12        17        15</Q>
+            <Q ID="10">        16        18        19        17</Q>
+            <Q ID="11">         5        21        22        20</Q>
+            <Q ID="12">        21        23        26         6</Q>
+            <Q ID="13">        22        24        27        23</Q>
+            <Q ID="14">        20         7        25        24</Q>
+            <Q ID="15">        11        26        27        25</Q>
+            <Q ID="16">        26        28        31        14</Q>
+            <Q ID="17">        27        29        32        28</Q>
+            <Q ID="18">        25        15        30        29</Q>
+            <Q ID="19">        19        31        32        30</Q>
+            <Q ID="20">        33        34        36         3</Q>
+            <Q ID="21">        33        35        38         4</Q>
+            <Q ID="22">        34        37        39        35</Q>
+            <Q ID="23">        36         6        40        37</Q>
+            <Q ID="24">        38        39        40        10</Q>
+            <Q ID="25">        38        41        43        13</Q>
+            <Q ID="26">        39        42        44        41</Q>
+            <Q ID="27">        40        14        45        42</Q>
+            <Q ID="28">        43        44        45        18</Q>
+            <Q ID="29">        36        46        47        21</Q>
+            <Q ID="30">        46        48        49        37</Q>
+            <Q ID="31">        47        23        50        48</Q>
+            <Q ID="32">        40        49        50        26</Q>
+            <Q ID="33">        49        51        52        42</Q>
+            <Q ID="34">        50        28        53        51</Q>
+            <Q ID="35">        45        52        53        31</Q>
+        </FACE>
+        <ELEMENT>
+            <H ID="0">    0     1     2     3     4     5 </H>
+            <H ID="1">    5     6     7     8     9    10 </H>
+            <H ID="2">   11     3    12    13    14    15 </H>
+            <H ID="3">   15     8    16    17    18    19 </H>
+            <H ID="4">   20    21    22    23     2    24 </H>
+            <H ID="5">   24    25    26    27     7    28 </H>
+            <H ID="6">   29    23    30    31    12    32 </H>
+            <H ID="7">   32    27    33    34    16    35 </H>
+        </ELEMENT>
+        <COMPOSITE>
+            <C ID="10"> H[0-2] </C>
+            <C ID="11"> H[3-4] </C>
+            <C ID="12"> H[5-6] </C>
+            <C ID="13"> H[7] </C>
+            <C ID="1"> F[0,11,20,29] </C>
+            <C ID="2"> F[1,6,21,25] </C>
+            <C ID="3"> F[22,26,30,33] </C>
+            <C ID="4"> F[13,17,31,34] </C>
+            <C ID="5"> F[4,9,14,18] </C>
+            <C ID="6"> F[10,19,28,35] </C>
+        </COMPOSITE>
+        <DOMAIN> C[10-13] </DOMAIN>
+    </GEOMETRY>
+    <EXPANSIONS>
+        <E COMPOSITE="C[10]" NUMMODES="4" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+        <E COMPOSITE="C[11]" NUMMODES="5" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+        <E COMPOSITE="C[12]" NUMMODES="6" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+        <E COMPOSITE="C[13]" NUMMODES="4" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+    </EXPANSIONS>
+    <CONDITIONS>
+        <SOLVERINFO>
+            <I PROPERTY="SolverType" VALUE="VelocityCorrectionScheme" />
+            <I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes" />
+            <I PROPERTY="AdvectionForm" VALUE="Convective" />
+            <I PROPERTY="Projection" VALUE="Galerkin" />
+            <I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder1" />
+        </SOLVERINFO>
+
+        <PARAMETERS>
+            <P> TimeStep      = 0.001 </P>
+            <P> NumSteps      = 10 </P>
+            <P> IO_CheckSteps = 10    </P>
+            <P> IO_InfoSteps  = 10    </P>
+            <P> Kinvis        = 1   </P>
+        </PARAMETERS>
+
+        <VARIABLES>
+            <V ID="0"> u </V>
+            <V ID="1"> v </V>
+            <V ID="2"> w </V>
+            <V ID="3"> p </V>
+        </VARIABLES>
+
+        <BOUNDARYREGIONS>
+            <B ID="0"> C[1] </B>  <!-- Inlet -->
+            <B ID="1"> C[6] </B>    <!-- Outlet -->
+            <B ID="2"> C[2] </B>    <!-- Wall -->
+            <B ID="3"> C[3] </B>    <!-- Wall left -->
+            <B ID="4"> C[4] </B>    <!-- Wall -->
+            <B ID="5"> C[5] </B>    <!-- Wall right -->
+        </BOUNDARYREGIONS>
+
+        <BOUNDARYCONDITIONS>
+            <REGION REF="0">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="1">
+                <N VAR="u" VALUE="0" />
+                <N VAR="v" VALUE="0" />
+                <N VAR="w" VALUE="0" />
+                <D VAR="p" VALUE="0" />
+            </REGION>
+            <REGION REF="2">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="3">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="4">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="5">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+        </BOUNDARYCONDITIONS>
+
+        <FUNCTION NAME="ExactSolution">
+            <E VAR="u" VALUE="0" />
+            <E VAR="v" VALUE="0" />
+            <E VAR="w" VALUE="y*(1-y)" />
+            <E VAR="p" VALUE="-2*Kinvis*(z-1)" />
+        </FUNCTION>
+        <FUNCTION NAME="InitialConditions">
+            <E VAR="u" VALUE="0" />
+            <E VAR="v" VALUE="0" />
+            <E VAR="w" VALUE="y*(1-y)" />
+            <E VAR="p" VALUE="-2*Kinvis*(z-1)" />
+        </FUNCTION>
+    </CONDITIONS>
+</NEKTAR>
diff --git a/solvers/IncNavierStokesSolver/Tests/Hex_channel_varP_par.tst b/solvers/IncNavierStokesSolver/Tests/Hex_channel_varP_par.tst
new file mode 100644
index 0000000000000000000000000000000000000000..4a60f5498737c6117e98d6cef874f53d958249e0
--- /dev/null
+++ b/solvers/IncNavierStokesSolver/Tests/Hex_channel_varP_par.tst
@@ -0,0 +1,24 @@
+<?xml version="1.0" encoding="utf-8"?>
+<test>
+    <description>3D channel flow, Hex elements, par(2), varP</description>
+    <executable>IncNavierStokesSolver</executable>
+    <parameters>-I GlobalSysSoln=XxtMultiLevelStaticCond Hex_channel_varP.xml</parameters>
+    <processes>2</processes>
+    <files>
+        <file description="Session File">Hex_channel_varP.xml</file>
+    </files>
+    <metrics>
+        <metric type="L2" id="1">
+            <value variable="u" tolerance="1e-12">0</value>
+            <value variable="v" tolerance="1e-12">0</value>
+            <value variable="w" tolerance="1e-12">0</value>
+            <value variable="p" tolerance="1e-11">0</value>
+        </metric>
+        <metric type="Linf" id="2">
+            <value variable="u" tolerance="1e-12">0</value>
+            <value variable="v" tolerance="1e-12">0</value>
+            <value variable="w" tolerance="1e-12">0</value>
+            <value variable="p" tolerance="1e-11">0</value>
+        </metric>
+    </metrics>
+</test>
diff --git a/solvers/IncNavierStokesSolver/Tests/Prism_channel_varP.tst b/solvers/IncNavierStokesSolver/Tests/Prism_channel_varP.tst
new file mode 100644
index 0000000000000000000000000000000000000000..928b5bb151ba60094b6a2fc362b5ad7e887df2cd
--- /dev/null
+++ b/solvers/IncNavierStokesSolver/Tests/Prism_channel_varP.tst
@@ -0,0 +1,23 @@
+<?xml version="1.0" encoding="utf-8"?>
+<test>
+    <description>3D channel flow, Prismatic elements, variable P</description>
+    <executable>IncNavierStokesSolver</executable>
+    <parameters>Prism_channel_varP.xml</parameters>
+    <files>
+        <file description="Session File">Prism_channel_varP.xml</file>
+    </files>
+    <metrics>
+        <metric type="L2" id="1">
+            <value variable="u" tolerance="1e-12">0</value>
+            <value variable="v" tolerance="1e-12">0</value>
+            <value variable="w" tolerance="1e-12">0</value>
+        <value variable="p" tolerance="1e-12">0</value>
+        </metric>
+        <metric type="Linf" id="2">
+            <value variable="u" tolerance="1e-12">0</value>
+            <value variable="v" tolerance="1e-12">0</value>
+            <value variable="w" tolerance="1e-12">0</value>
+        <value variable="p" tolerance="1e-12">0</value>
+        </metric>
+    </metrics>
+</test>
diff --git a/solvers/IncNavierStokesSolver/Tests/Prism_channel_varP.xml b/solvers/IncNavierStokesSolver/Tests/Prism_channel_varP.xml
new file mode 100644
index 0000000000000000000000000000000000000000..05a2df29be074dec13c6674aa85719790d8ef19e
--- /dev/null
+++ b/solvers/IncNavierStokesSolver/Tests/Prism_channel_varP.xml
@@ -0,0 +1,280 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <GEOMETRY DIM="3" SPACE="3">
+        <VERTEX>
+            <V ID="0">0.00000000e+00 0.00000000e+00 5.00000000e-01</V>
+            <V ID="1">5.00000000e-01 0.00000000e+00 5.00000000e-01</V>
+            <V ID="2">5.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="3">0.00000000e+00 0.00000000e+00 0.00000000e+00</V>
+            <V ID="4">5.00000000e-01 5.00000000e-01 5.00000000e-01</V>
+            <V ID="5">5.00000000e-01 5.00000000e-01 0.00000000e+00</V>
+            <V ID="6">0.00000000e+00 0.00000000e+00 1.00000000e+00</V>
+            <V ID="7">5.00000000e-01 0.00000000e+00 1.00000000e+00</V>
+            <V ID="8">5.00000000e-01 5.00000000e-01 1.00000000e+00</V>
+            <V ID="9">0.00000000e+00 5.00000000e-01 5.00000000e-01</V>
+            <V ID="10">0.00000000e+00 5.00000000e-01 0.00000000e+00</V>
+            <V ID="11">0.00000000e+00 5.00000000e-01 1.00000000e+00</V>
+            <V ID="12">5.00000000e-01 1.00000000e+00 5.00000000e-01</V>
+            <V ID="13">5.00000000e-01 1.00000000e+00 0.00000000e+00</V>
+            <V ID="14">5.00000000e-01 1.00000000e+00 1.00000000e+00</V>
+            <V ID="15">0.00000000e+00 1.00000000e+00 5.00000000e-01</V>
+            <V ID="16">0.00000000e+00 1.00000000e+00 0.00000000e+00</V>
+            <V ID="17">0.00000000e+00 1.00000000e+00 1.00000000e+00</V>
+            <V ID="18">1.00000000e+00 0.00000000e+00 5.00000000e-01</V>
+            <V ID="19">1.00000000e+00 0.00000000e+00 0.00000000e+00</V>
+            <V ID="20">1.00000000e+00 5.00000000e-01 5.00000000e-01</V>
+            <V ID="21">1.00000000e+00 5.00000000e-01 0.00000000e+00</V>
+            <V ID="22">1.00000000e+00 0.00000000e+00 1.00000000e+00</V>
+            <V ID="23">1.00000000e+00 5.00000000e-01 1.00000000e+00</V>
+            <V ID="24">1.00000000e+00 1.00000000e+00 5.00000000e-01</V>
+            <V ID="25">1.00000000e+00 1.00000000e+00 0.00000000e+00</V>
+            <V ID="26">1.00000000e+00 1.00000000e+00 1.00000000e+00</V>
+        </VERTEX>
+        <EDGE>
+            <E ID="0">    0  1   </E>
+            <E ID="1">    0  3   </E>
+            <E ID="2">    0  4   </E>
+            <E ID="3">    1  2   </E>
+            <E ID="4">    1  4   </E>
+            <E ID="5">    2  5   </E>
+            <E ID="6">    3  2   </E>
+            <E ID="7">    3  5   </E>
+            <E ID="8">    4  5   </E>
+            <E ID="9">    6  7   </E>
+            <E ID="10">    6  0   </E>
+            <E ID="11">    6  8   </E>
+            <E ID="12">    7  1   </E>
+            <E ID="13">    7  8   </E>
+            <E ID="14">    8  4   </E>
+            <E ID="15">    0  9   </E>
+            <E ID="16">    4  9   </E>
+            <E ID="17">    5  10   </E>
+            <E ID="18">    3  10   </E>
+            <E ID="19">    9  10   </E>
+            <E ID="20">    6  11   </E>
+            <E ID="21">    8  11   </E>
+            <E ID="22">   11  9   </E>
+            <E ID="23">    9  12   </E>
+            <E ID="24">    4  12   </E>
+            <E ID="25">    5  13   </E>
+            <E ID="26">   10  13   </E>
+            <E ID="27">   12  13   </E>
+            <E ID="28">   11  14   </E>
+            <E ID="29">    8  14   </E>
+            <E ID="30">   14  12   </E>
+            <E ID="31">    9  15   </E>
+            <E ID="32">   12  15   </E>
+            <E ID="33">   13  16   </E>
+            <E ID="34">   10  16   </E>
+            <E ID="35">   15  16   </E>
+            <E ID="36">   11  17   </E>
+            <E ID="37">   14  17   </E>
+            <E ID="38">   17  15   </E>
+            <E ID="39">    1  18   </E>
+            <E ID="40">    1  20   </E>
+            <E ID="41">   18  19   </E>
+            <E ID="42">   18  20   </E>
+            <E ID="43">   19  21   </E>
+            <E ID="44">    2  19   </E>
+            <E ID="45">    2  21   </E>
+            <E ID="46">   20  21   </E>
+            <E ID="47">    7  22   </E>
+            <E ID="48">    7  23   </E>
+            <E ID="49">   22  18   </E>
+            <E ID="50">   22  23   </E>
+            <E ID="51">   23  20   </E>
+            <E ID="52">   20  4   </E>
+            <E ID="53">   21  5   </E>
+            <E ID="54">   23  8   </E>
+            <E ID="55">    4  24   </E>
+            <E ID="56">   20  24   </E>
+            <E ID="57">   21  25   </E>
+            <E ID="58">    5  25   </E>
+            <E ID="59">   24  25   </E>
+            <E ID="60">    8  26   </E>
+            <E ID="61">   23  26   </E>
+            <E ID="62">   26  24   </E>
+            <E ID="63">   24  12   </E>
+            <E ID="64">   25  13   </E>
+            <E ID="65">   26  14   </E>
+        </EDGE>
+        <FACE>
+            <Q ID="0">         0         3         6         1</Q>
+            <T ID="1">         0         4         2</T>
+            <Q ID="2">         3         5         8         4</Q>
+            <T ID="3">         6         5         7</T>
+            <Q ID="4">         1         7         8         2</Q>
+            <Q ID="5">         9        12         0        10</Q>
+            <T ID="6">         9        13        11</T>
+            <Q ID="7">        12         4        14        13</Q>
+            <Q ID="8">        10         2        14        11</Q>
+            <Q ID="9">        15         1        18        19</Q>
+            <T ID="10">        15         2        16</T>
+            <T ID="11">        18         7        17</T>
+            <Q ID="12">        19        17         8        16</Q>
+            <Q ID="13">        20        10        15        22</Q>
+            <T ID="14">        20        11        21</T>
+            <Q ID="15">        22        16        14        21</Q>
+            <Q ID="16">        23        19        26        27</Q>
+            <T ID="17">        23        16        24</T>
+            <T ID="18">        26        17        25</T>
+            <Q ID="19">        27        25         8        24</Q>
+            <Q ID="20">        28        22        23        30</Q>
+            <T ID="21">        28        21        29</T>
+            <Q ID="22">        30        24        14        29</Q>
+            <Q ID="23">        32        35        33        27</Q>
+            <T ID="24">        32        31        23</T>
+            <Q ID="25">        35        34        19        31</Q>
+            <T ID="26">        33        34        26</T>
+            <Q ID="27">        37        38        32        30</Q>
+            <T ID="28">        37        36        28</T>
+            <Q ID="29">        38        31        22        36</Q>
+            <Q ID="30">        39        41        44         3</Q>
+            <T ID="31">        39        42        40</T>
+            <Q ID="32">        41        43        46        42</Q>
+            <T ID="33">        44        43        45</T>
+            <Q ID="34">         3        45        46        40</Q>
+            <Q ID="35">        47        49        39        12</Q>
+            <T ID="36">        47        50        48</T>
+            <Q ID="37">        49        42        51        50</Q>
+            <Q ID="38">        12        40        51        48</Q>
+            <T ID="39">        40        52         4</T>
+            <Q ID="40">        46        53         8        52</Q>
+            <T ID="41">        45        53         5</T>
+            <T ID="42">        48        54        13</T>
+            <Q ID="43">        51        52        14        54</Q>
+            <Q ID="44">        56        59        57        46</Q>
+            <T ID="45">        56        55        52</T>
+            <Q ID="46">        59        58         8        55</Q>
+            <T ID="47">        57        58        53</T>
+            <Q ID="48">        61        62        56        51</Q>
+            <T ID="49">        61        60        54</T>
+            <Q ID="50">        62        55        14        60</Q>
+            <Q ID="51">        63        27        64        59</Q>
+            <T ID="52">        63        24        55</T>
+            <T ID="53">        64        25        58</T>
+            <Q ID="54">        65        30        63        62</Q>
+            <T ID="55">        65        29        60</T>
+        </FACE>
+        <ELEMENT>
+            <R ID="0">    0     1     2     3     4 </R>
+            <R ID="1">    5     6     7     1     8 </R>
+            <R ID="2">    9    10     4    11    12 </R>
+            <R ID="3">   13    14     8    10    15 </R>
+            <R ID="4">   16    17    12    18    19 </R>
+            <R ID="5">   20    21    15    17    22 </R>
+            <R ID="6">   23    24    25    26    16 </R>
+            <R ID="7">   27    28    29    24    20 </R>
+            <R ID="8">   30    31    32    33    34 </R>
+            <R ID="9">   35    36    37    31    38 </R>
+            <R ID="10">   34    39    40    41     2 </R>
+            <R ID="11">   38    42    43    39     7 </R>
+            <R ID="12">   44    45    46    47    40 </R>
+            <R ID="13">   48    49    50    45    43 </R>
+            <R ID="14">   51    52    19    53    46 </R>
+            <R ID="15">   54    55    22    52    50 </R>
+        </ELEMENT>
+        <COMPOSITE>
+            <C ID="10"> R[0-5] </C>
+            <C ID="11"> R[6-10] </C>
+            <C ID="12"> R[11-15] </C>
+            <C ID="1"> F[3,11,18,26,33,41,47,53] </C>
+            <C ID="2"> F[0,5,30,35] </C>
+            <C ID="3"> F[32,37,44,48] </C>
+            <C ID="4"> F[23,27,51,54] </C>
+            <C ID="5"> F[9,13,25,29] </C>
+            <C ID="6"> F[6,14,21,28,36,42,49,55] </C>
+        </COMPOSITE>
+        <DOMAIN> C[10-12] </DOMAIN>
+    </GEOMETRY>
+    <EXPANSIONS>
+        <E COMPOSITE="C[10]" NUMMODES="4" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+        <E COMPOSITE="C[11]" NUMMODES="6" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+        <E COMPOSITE="C[12]" NUMMODES="5" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+    </EXPANSIONS>
+    <CONDITIONS>
+        <SOLVERINFO>
+            <I PROPERTY="SolverType" VALUE="VelocityCorrectionScheme" />
+            <I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes" />
+            <I PROPERTY="AdvectionForm" VALUE="Convective" />
+            <I PROPERTY="Projection" VALUE="Galerkin" />
+            <I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder1" />
+        </SOLVERINFO>
+
+        <PARAMETERS>
+            <P> TimeStep      = 0.001 </P>
+            <P> NumSteps      = 10 </P>
+            <P> IO_CheckSteps = 10    </P>
+            <P> IO_InfoSteps  = 10    </P>
+            <P> Kinvis        = 1   </P>
+        </PARAMETERS>
+
+        <VARIABLES>
+            <V ID="0"> u </V>
+            <V ID="1"> v </V>
+            <V ID="2"> w </V>
+            <V ID="3"> p </V>
+        </VARIABLES>
+
+        <BOUNDARYREGIONS>
+            <B ID="0"> C[1] </B>  <!-- Inlet -->
+            <B ID="1"> C[6] </B>    <!-- Outlet -->
+            <B ID="2"> C[2] </B>    <!-- Wall -->
+            <B ID="3"> C[3] </B>    <!-- Wall left -->
+            <B ID="4"> C[4] </B>    <!-- Wall -->
+            <B ID="5"> C[5] </B>    <!-- Wall right -->
+        </BOUNDARYREGIONS>
+
+        <BOUNDARYCONDITIONS>
+            <REGION REF="0">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="1">
+                <N VAR="u" VALUE="0" />
+                <N VAR="v" VALUE="0" />
+                <N VAR="w" VALUE="0" />
+                <D VAR="p" VALUE="0" />
+            </REGION>
+            <REGION REF="2">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="3">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="4">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="5">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+        </BOUNDARYCONDITIONS>
+
+        <FUNCTION NAME="ExactSolution">
+            <E VAR="u" VALUE="0" />
+            <E VAR="v" VALUE="0" />
+            <E VAR="w" VALUE="y*(1-y)" />
+            <E VAR="p" VALUE="-2*Kinvis*(z-1)" />
+        </FUNCTION>
+        <FUNCTION NAME="InitialConditions">
+            <E VAR="u" VALUE="0" />
+            <E VAR="v" VALUE="0" />
+            <E VAR="w" VALUE="y*(1-y)" />
+            <E VAR="p" VALUE="-2*Kinvis*(z-1)" />
+        </FUNCTION>
+    </CONDITIONS>
+</NEKTAR>
diff --git a/solvers/IncNavierStokesSolver/Tests/Pyr_channel_varP.tst b/solvers/IncNavierStokesSolver/Tests/Pyr_channel_varP.tst
new file mode 100644
index 0000000000000000000000000000000000000000..4cdc04b9d8afc900d20b1dabfec1aabdfa6b7b48
--- /dev/null
+++ b/solvers/IncNavierStokesSolver/Tests/Pyr_channel_varP.tst
@@ -0,0 +1,23 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<test>
+    <description>3D channel flow, Pyramidic elements, variable P</description>
+    <executable>IncNavierStokesSolver</executable>
+    <parameters>Pyr_channel_varP.xml</parameters>
+    <files>
+        <file description="Session File">Pyr_channel_varP.xml</file>
+    </files>
+    <metrics>
+        <metric type="L2" id="1">
+            <value variable="u" tolerance="1e-12">0</value>
+            <value variable="v" tolerance="1e-12">0</value>
+            <value variable="w" tolerance="1e-12">0</value>
+        <value variable="p" tolerance="1e-12">0</value>
+        </metric>
+        <metric type="Linf" id="2">
+            <value variable="u" tolerance="1e-12">0</value>
+            <value variable="v" tolerance="1e-12">0</value>
+            <value variable="w" tolerance="1e-12">0</value>
+        <value variable="p" tolerance="1e-12">0</value>
+        </metric>
+    </metrics>
+</test>
diff --git a/solvers/IncNavierStokesSolver/Tests/Pyr_channel_varP.xml b/solvers/IncNavierStokesSolver/Tests/Pyr_channel_varP.xml
new file mode 100644
index 0000000000000000000000000000000000000000..25b08d337fb096b0d8dc5ee3781ca9725f3a609c
--- /dev/null
+++ b/solvers/IncNavierStokesSolver/Tests/Pyr_channel_varP.xml
@@ -0,0 +1,168 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <EXPANSIONS>
+        <E COMPOSITE="C[10]" NUMMODES="4" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+        <E COMPOSITE="C[11]" NUMMODES="6" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+        <E COMPOSITE="C[12]" NUMMODES="5" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+    </EXPANSIONS>
+    <CONDITIONS>
+        <SOLVERINFO>
+            <I PROPERTY="SolverType" VALUE="VelocityCorrectionScheme" />
+            <I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes" />
+            <I PROPERTY="AdvectionForm" VALUE="Convective" />
+            <I PROPERTY="Projection" VALUE="Galerkin" />
+            <I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder1" />
+        </SOLVERINFO>
+
+        <PARAMETERS>
+            <P> TimeStep      = 0.001 </P>
+            <P> NumSteps      = 100 </P>
+            <P> IO_InfoSteps  = 100    </P>
+            <P> Kinvis        = 1   </P>
+        </PARAMETERS>
+
+        <VARIABLES>
+            <V ID="0"> u </V>
+            <V ID="1"> v </V>
+            <V ID="2"> w </V>
+            <V ID="3"> p </V>
+        </VARIABLES>
+
+        <BOUNDARYREGIONS>
+            <B ID="0"> C[1] </B>  <!-- Inlet -->
+            <B ID="1"> C[6] </B>    <!-- Outlet -->
+            <B ID="2"> C[2] </B>    <!-- Wall -->
+            <B ID="3"> C[3] </B>    <!-- Wall left -->
+            <B ID="4"> C[4] </B>    <!-- Wall -->
+            <B ID="5"> C[5] </B>    <!-- Wall right -->
+        </BOUNDARYREGIONS>
+
+        <BOUNDARYCONDITIONS>
+            <REGION REF="0">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="1">
+                <N VAR="u" VALUE="0" />
+                <N VAR="v" VALUE="0" />
+                <N VAR="w" VALUE="0" />
+                <D VAR="p" VALUE="0" />
+            </REGION>
+            <REGION REF="2">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="3">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="4">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="5">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+        </BOUNDARYCONDITIONS>
+
+        <FUNCTION NAME="ExactSolution">
+            <E VAR="u" VALUE="0" />
+            <E VAR="v" VALUE="0" />
+            <E VAR="w" VALUE="y*(1-y)" />
+            <E VAR="p" VALUE="-2*Kinvis*(z-1)" />
+        </FUNCTION>
+
+        <FUNCTION NAME="InitialConditions">
+            <E VAR="u" VALUE="0" />
+            <E VAR="v" VALUE="0" />
+            <E VAR="w" VALUE="y*(1-y)" />
+            <E VAR="p" VALUE="-2*Kinvis*(z-1)" />
+        </FUNCTION>
+    </CONDITIONS>
+    <GEOMETRY DIM="3" SPACE="3">
+        <VERTEX>
+            <V ID="0">0.00000000e+00 0.00000000e+00 0.00000000e+00</V>
+            <V ID="1">1.00000000e+00 0.00000000e+00 0.00000000e+00</V>
+            <V ID="2">1.00000000e+00 1.00000000e+00 0.00000000e+00</V>
+            <V ID="3">0.00000000e+00 1.00000000e+00 0.00000000e+00</V>
+            <V ID="4">5.00000000e-01 5.00000000e-01 5.00000000e-01</V>
+            <V ID="5">0.00000000e+00 0.00000000e+00 1.00000000e+00</V>
+            <V ID="6">1.00000000e+00 0.00000000e+00 1.00000000e+00</V>
+            <V ID="7">1.00000000e+00 1.00000000e+00 1.00000000e+00</V>
+            <V ID="8">0.00000000e+00 1.00000000e+00 1.00000000e+00</V>
+        </VERTEX>
+        <EDGE>
+            <E ID="0">    0  1   </E>
+            <E ID="1">    1  2   </E>
+            <E ID="2">    3  2   </E>
+            <E ID="3">    0  3   </E>
+            <E ID="4">    0  4   </E>
+            <E ID="5">    1  4   </E>
+            <E ID="6">    2  4   </E>
+            <E ID="7">    3  4   </E>
+            <E ID="8">    0  5   </E>
+            <E ID="9">    5  6   </E>
+            <E ID="10">    1  6   </E>
+            <E ID="11">    5  4   </E>
+            <E ID="12">    6  4   </E>
+            <E ID="13">    6  7   </E>
+            <E ID="14">    2  7   </E>
+            <E ID="15">    7  4   </E>
+            <E ID="16">    7  8   </E>
+            <E ID="17">    3  8   </E>
+            <E ID="18">    8  4   </E>
+            <E ID="19">    8  5   </E>
+        </EDGE>
+        <FACE>
+            <Q ID="0">         0         1         2         3</Q>
+            <T ID="1">         0         5         4</T>
+            <T ID="2">         1         6         5</T>
+            <T ID="3">         2         6         7</T>
+            <T ID="4">         3         7         4</T>
+            <Q ID="5">         8         9        10         0</Q>
+            <T ID="6">         8        11         4</T>
+            <T ID="7">         9        12        11</T>
+            <T ID="8">        10        12         5</T>
+            <Q ID="9">        10        13        14         1</Q>
+            <T ID="10">        13        15        12</T>
+            <T ID="11">        14        15         6</T>
+            <Q ID="12">        14        16        17         2</Q>
+            <T ID="13">        16        18        15</T>
+            <T ID="14">        17        18         7</T>
+            <Q ID="15">        17        19         8         3</Q>
+            <T ID="16">        19        11        18</T>
+            <Q ID="17">        19        16        13         9</Q>
+        </FACE>
+        <ELEMENT>
+            <P ID="0">    0     1     2     3     4 </P>
+            <P ID="1">    5     6     7     8     1 </P>
+            <P ID="2">    9     8    10    11     2 </P>
+            <P ID="3">   12    11    13    14     3 </P>
+            <P ID="4">   15    14    16     6     4 </P>
+            <P ID="5">   17    16    13    10     7 </P>
+        </ELEMENT>
+        <COMPOSITE>
+            <C ID="10"> P[0-2] </C>
+            <C ID="11"> P[3-4] </C>
+            <C ID="12"> P[5] </C>
+            <C ID="1"> F[0] </C>
+            <C ID="2"> F[5] </C>
+            <C ID="3"> F[9] </C>
+            <C ID="4"> F[12] </C>
+            <C ID="5"> F[15] </C>
+            <C ID="6"> F[17] </C>
+        </COMPOSITE>
+        <DOMAIN> C[10-12] </DOMAIN>
+    </GEOMETRY>
+</NEKTAR>
diff --git a/solvers/IncNavierStokesSolver/Tests/Tet_channel_varP.tst b/solvers/IncNavierStokesSolver/Tests/Tet_channel_varP.tst
new file mode 100644
index 0000000000000000000000000000000000000000..edb538e003081216d889dc4283bd64b94ddbf6d2
--- /dev/null
+++ b/solvers/IncNavierStokesSolver/Tests/Tet_channel_varP.tst
@@ -0,0 +1,24 @@
+<?xml version="1.0" encoding="utf-8"?>
+<test>
+    <description>3D channel flow, Tetrahedral elements, variable P</description>
+    <executable>IncNavierStokesSolver</executable>
+    <parameters>Tet_channel_varP.xml</parameters>
+    <files>
+        <file description="Session File">Tet_channel_varP.xml</file>
+    </files>
+    <metrics>
+        <metric type="L2" id="1">
+            <value variable="u" tolerance="1e-12">0</value>
+            <value variable="v" tolerance="1e-12">0</value>
+            <value variable="w" tolerance="1e-12">0</value>
+        <value variable="p" tolerance="1e-12">0</value>
+        </metric>
+        <metric type="Linf" id="2">
+            <value variable="u" tolerance="1e-12">0</value>
+            <value variable="v" tolerance="1e-12">0</value>
+            <value variable="w" tolerance="1e-12">0</value>
+        <value variable="p" tolerance="1e-12">0</value>
+        </metric>
+    </metrics>
+</test>
+
diff --git a/solvers/IncNavierStokesSolver/Tests/Tet_channel_varP.xml b/solvers/IncNavierStokesSolver/Tests/Tet_channel_varP.xml
new file mode 100644
index 0000000000000000000000000000000000000000..c3d55e4f44d38bdb3e09daba6978c817ddaa1202
--- /dev/null
+++ b/solvers/IncNavierStokesSolver/Tests/Tet_channel_varP.xml
@@ -0,0 +1,411 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <GEOMETRY DIM="3" SPACE="3">
+        <VERTEX>
+            <V ID="0">0.00000000e+00 0.00000000e+00 0.00000000e+00</V>
+            <V ID="1">5.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="2">5.00000000e-01 5.00000000e-01 0.00000000e+00</V>
+            <V ID="3">5.00000000e-01 0.00000000e+00 5.00000000e-01</V>
+            <V ID="4">0.00000000e+00 0.00000000e+00 5.00000000e-01</V>
+            <V ID="5">5.00000000e-01 5.00000000e-01 5.00000000e-01</V>
+            <V ID="6">0.00000000e+00 0.00000000e+00 1.00000000e+00</V>
+            <V ID="7">5.00000000e-01 0.00000000e+00 1.00000000e+00</V>
+            <V ID="8">5.00000000e-01 5.00000000e-01 1.00000000e+00</V>
+            <V ID="9">0.00000000e+00 5.00000000e-01 0.00000000e+00</V>
+            <V ID="10">0.00000000e+00 5.00000000e-01 5.00000000e-01</V>
+            <V ID="11">0.00000000e+00 5.00000000e-01 1.00000000e+00</V>
+            <V ID="12">5.00000000e-01 1.00000000e+00 0.00000000e+00</V>
+            <V ID="13">5.00000000e-01 1.00000000e+00 5.00000000e-01</V>
+            <V ID="14">5.00000000e-01 1.00000000e+00 1.00000000e+00</V>
+            <V ID="15">0.00000000e+00 1.00000000e+00 0.00000000e+00</V>
+            <V ID="16">0.00000000e+00 1.00000000e+00 5.00000000e-01</V>
+            <V ID="17">0.00000000e+00 1.00000000e+00 1.00000000e+00</V>
+            <V ID="18">1.00000000e+00 0.00000000e+00 0.00000000e+00</V>
+            <V ID="19">1.00000000e+00 5.00000000e-01 0.00000000e+00</V>
+            <V ID="20">1.00000000e+00 0.00000000e+00 5.00000000e-01</V>
+            <V ID="21">1.00000000e+00 5.00000000e-01 5.00000000e-01</V>
+            <V ID="22">1.00000000e+00 0.00000000e+00 1.00000000e+00</V>
+            <V ID="23">1.00000000e+00 5.00000000e-01 1.00000000e+00</V>
+            <V ID="24">1.00000000e+00 1.00000000e+00 0.00000000e+00</V>
+            <V ID="25">1.00000000e+00 1.00000000e+00 5.00000000e-01</V>
+            <V ID="26">1.00000000e+00 1.00000000e+00 1.00000000e+00</V>
+        </VERTEX>
+        <EDGE>
+            <E ID="0">    0  1   </E>
+            <E ID="1">    0  2   </E>
+            <E ID="2">    0  3   </E>
+            <E ID="3">    1  2   </E>
+            <E ID="4">    1  3   </E>
+            <E ID="5">    2  3   </E>
+            <E ID="6">    3  4   </E>
+            <E ID="7">    3  5   </E>
+            <E ID="8">    4  5   </E>
+            <E ID="9">    4  0   </E>
+            <E ID="10">    5  0   </E>
+            <E ID="11">    2  5   </E>
+            <E ID="12">    4  6   </E>
+            <E ID="13">    3  6   </E>
+            <E ID="14">    5  6   </E>
+            <E ID="15">    7  6   </E>
+            <E ID="16">    7  8   </E>
+            <E ID="17">    7  3   </E>
+            <E ID="18">    6  8   </E>
+            <E ID="19">    8  3   </E>
+            <E ID="20">    8  5   </E>
+            <E ID="21">    0  9   </E>
+            <E ID="22">    0  10   </E>
+            <E ID="23">    2  9   </E>
+            <E ID="24">    2  10   </E>
+            <E ID="25">    9  10   </E>
+            <E ID="26">    5  10   </E>
+            <E ID="27">    4  10   </E>
+            <E ID="28">   10  6   </E>
+            <E ID="29">    8  11   </E>
+            <E ID="30">    8  10   </E>
+            <E ID="31">    6  11   </E>
+            <E ID="32">   11  10   </E>
+            <E ID="33">    9  12   </E>
+            <E ID="34">    9  13   </E>
+            <E ID="35">    2  12   </E>
+            <E ID="36">    2  13   </E>
+            <E ID="37">   12  13   </E>
+            <E ID="38">    5  13   </E>
+            <E ID="39">   10  13   </E>
+            <E ID="40">   13  8   </E>
+            <E ID="41">    8  14   </E>
+            <E ID="42">   11  14   </E>
+            <E ID="43">   11  13   </E>
+            <E ID="44">   14  13   </E>
+            <E ID="45">    9  15   </E>
+            <E ID="46">   12  15   </E>
+            <E ID="47">   15  13   </E>
+            <E ID="48">   13  16   </E>
+            <E ID="49">   10  16   </E>
+            <E ID="50">   10  15   </E>
+            <E ID="51">   16  15   </E>
+            <E ID="52">   10  17   </E>
+            <E ID="53">   13  17   </E>
+            <E ID="54">   16  17   </E>
+            <E ID="55">   14  17   </E>
+            <E ID="56">   11  17   </E>
+            <E ID="57">    1  18   </E>
+            <E ID="58">    1  19   </E>
+            <E ID="59">   18  19   </E>
+            <E ID="60">   18  3   </E>
+            <E ID="61">   19  3   </E>
+            <E ID="62">   20  3   </E>
+            <E ID="63">   20  21   </E>
+            <E ID="64">   20  18   </E>
+            <E ID="65">    3  21   </E>
+            <E ID="66">   21  18   </E>
+            <E ID="67">   21  19   </E>
+            <E ID="68">    3  22   </E>
+            <E ID="69">   20  22   </E>
+            <E ID="70">   21  22   </E>
+            <E ID="71">   22  7   </E>
+            <E ID="72">   22  23   </E>
+            <E ID="73">    7  23   </E>
+            <E ID="74">   23  3   </E>
+            <E ID="75">   21  23   </E>
+            <E ID="76">   19  2   </E>
+            <E ID="77">   21  5   </E>
+            <E ID="78">   21  2   </E>
+            <E ID="79">   21  8   </E>
+            <E ID="80">   23  8   </E>
+            <E ID="81">    2  24   </E>
+            <E ID="82">   19  24   </E>
+            <E ID="83">   24  21   </E>
+            <E ID="84">   21  25   </E>
+            <E ID="85">    5  25   </E>
+            <E ID="86">    5  24   </E>
+            <E ID="87">   25  24   </E>
+            <E ID="88">    5  26   </E>
+            <E ID="89">   21  26   </E>
+            <E ID="90">   25  26   </E>
+            <E ID="91">   23  26   </E>
+            <E ID="92">    8  26   </E>
+            <E ID="93">   24  12   </E>
+            <E ID="94">   24  13   </E>
+            <E ID="95">   25  13   </E>
+            <E ID="96">   13  26   </E>
+            <E ID="97">   26  14   </E>
+        </EDGE>
+        <FACE>
+            <T ID="0">         0         3         1</T>
+            <T ID="1">         0         4         2</T>
+            <T ID="2">         3         5         4</T>
+            <T ID="3">         1         5         2</T>
+            <T ID="4">         9         6         2</T>
+            <T ID="5">         9         8        10</T>
+            <T ID="6">         6         7         8</T>
+            <T ID="7">         2         7        10</T>
+            <T ID="8">         1        10        11</T>
+            <T ID="9">         5         7        11</T>
+            <T ID="10">        12        13         6</T>
+            <T ID="11">        12        14         8</T>
+            <T ID="12">        13         7        14</T>
+            <T ID="13">        15        17        13</T>
+            <T ID="14">        15        16        18</T>
+            <T ID="15">        17        19        16</T>
+            <T ID="16">        13        19        18</T>
+            <T ID="17">        18        20        14</T>
+            <T ID="18">        19        20         7</T>
+            <T ID="19">        21         1        23</T>
+            <T ID="20">        21        22        25</T>
+            <T ID="21">         1        24        22</T>
+            <T ID="22">        23        24        25</T>
+            <T ID="23">         9        22        27</T>
+            <T ID="24">        22        26        10</T>
+            <T ID="25">        27        26         8</T>
+            <T ID="26">        24        26        11</T>
+            <T ID="27">        12        27        28</T>
+            <T ID="28">        28        26        14</T>
+            <T ID="29">        31        28        32</T>
+            <T ID="30">        31        18        29</T>
+            <T ID="31">        28        30        18</T>
+            <T ID="32">        32        30        29</T>
+            <T ID="33">        30        20        26</T>
+            <T ID="34">        33        23        35</T>
+            <T ID="35">        33        34        37</T>
+            <T ID="36">        23        36        34</T>
+            <T ID="37">        35        36        37</T>
+            <T ID="38">        36        39        24</T>
+            <T ID="39">        36        38        11</T>
+            <T ID="40">        39        26        38</T>
+            <T ID="41">        34        39        25</T>
+            <T ID="42">        39        40        30</T>
+            <T ID="43">        40        20        38</T>
+            <T ID="44">        42        43        44</T>
+            <T ID="45">        42        29        41</T>
+            <T ID="46">        43        40        29</T>
+            <T ID="47">        44        40        41</T>
+            <T ID="48">        43        32        39</T>
+            <T ID="49">        46        45        33</T>
+            <T ID="50">        46        47        37</T>
+            <T ID="51">        45        34        47</T>
+            <T ID="52">        51        47        48</T>
+            <T ID="53">        51        50        49</T>
+            <T ID="54">        47        39        50</T>
+            <T ID="55">        48        39        49</T>
+            <T ID="56">        45        25        50</T>
+            <T ID="57">        54        48        53</T>
+            <T ID="58">        54        49        52</T>
+            <T ID="59">        53        39        52</T>
+            <T ID="60">        55        42        56</T>
+            <T ID="61">        55        44        53</T>
+            <T ID="62">        56        43        53</T>
+            <T ID="63">        56        52        32</T>
+            <T ID="64">        57        59        58</T>
+            <T ID="65">        57        60         4</T>
+            <T ID="66">        59        61        60</T>
+            <T ID="67">        58        61         4</T>
+            <T ID="68">        64        60        62</T>
+            <T ID="69">        64        66        63</T>
+            <T ID="70">        60        65        66</T>
+            <T ID="71">        62        65        63</T>
+            <T ID="72">        59        67        66</T>
+            <T ID="73">        61        65        67</T>
+            <T ID="74">        69        62        68</T>
+            <T ID="75">        69        63        70</T>
+            <T ID="76">        68        65        70</T>
+            <T ID="77">        71        73        72</T>
+            <T ID="78">        71        17        68</T>
+            <T ID="79">        73        74        17</T>
+            <T ID="80">        72        74        68</T>
+            <T ID="81">        72        70        75</T>
+            <T ID="82">        74        65        75</T>
+            <T ID="83">        58        76         3</T>
+            <T ID="84">        76         5        61</T>
+            <T ID="85">         5        65        78</T>
+            <T ID="86">        65        77         7</T>
+            <T ID="87">        78        77        11</T>
+            <T ID="88">        76        78        67</T>
+            <T ID="89">        65        19        79</T>
+            <T ID="90">        79        20        77</T>
+            <T ID="91">        73        80        16</T>
+            <T ID="92">        74        19        80</T>
+            <T ID="93">        75        79        80</T>
+            <T ID="94">        82        81        76</T>
+            <T ID="95">        82        83        67</T>
+            <T ID="96">        81        78        83</T>
+            <T ID="97">        87        83        84</T>
+            <T ID="98">        87        86        85</T>
+            <T ID="99">        83        77        86</T>
+            <T ID="100">        84        77        85</T>
+            <T ID="101">        81        11        86</T>
+            <T ID="102">        90        84        89</T>
+            <T ID="103">        90        85        88</T>
+            <T ID="104">        89        77        88</T>
+            <T ID="105">        91        89        75</T>
+            <T ID="106">        91        92        80</T>
+            <T ID="107">        89        79        92</T>
+            <T ID="108">        92        20        88</T>
+            <T ID="109">        93        35        81</T>
+            <T ID="110">        93        37        94</T>
+            <T ID="111">        81        36        94</T>
+            <T ID="112">        87        95        94</T>
+            <T ID="113">        95        38        85</T>
+            <T ID="114">        94        38        86</T>
+            <T ID="115">        90        96        95</T>
+            <T ID="116">        96        38        88</T>
+            <T ID="117">        97        44        96</T>
+            <T ID="118">        97        41        92</T>
+            <T ID="119">        96        40        92</T>
+        </FACE>
+        <ELEMENT>
+            <A ID="0">    0     1     2     3 </A>
+            <A ID="1">    4     5     6     7 </A>
+            <A ID="2">    3     8     7     9 </A>
+            <A ID="3">   10    11    12     6 </A>
+            <A ID="4">   13    14    15    16 </A>
+            <A ID="5">   16    12    17    18 </A>
+            <A ID="6">   19    20    21    22 </A>
+            <A ID="7">   23     5    24    25 </A>
+            <A ID="8">   21     8    26    24 </A>
+            <A ID="9">   27    11    25    28 </A>
+            <A ID="10">   29    30    31    32 </A>
+            <A ID="11">   31    28    33    17 </A>
+            <A ID="12">   34    35    36    37 </A>
+            <A ID="13">   38    39    40    26 </A>
+            <A ID="14">   36    22    38    41 </A>
+            <A ID="15">   42    40    43    33 </A>
+            <A ID="16">   44    45    46    47 </A>
+            <A ID="17">   48    46    32    42 </A>
+            <A ID="18">   49    50    51    35 </A>
+            <A ID="19">   52    53    54    55 </A>
+            <A ID="20">   51    56    41    54 </A>
+            <A ID="21">   57    58    55    59 </A>
+            <A ID="22">   60    61    44    62 </A>
+            <A ID="23">   62    63    59    48 </A>
+            <A ID="24">   64    65    66    67 </A>
+            <A ID="25">   68    69    70    71 </A>
+            <A ID="26">   66    72    73    70 </A>
+            <A ID="27">   74    75    71    76 </A>
+            <A ID="28">   77    78    79    80 </A>
+            <A ID="29">   80    81    76    82 </A>
+            <A ID="30">   83    67    84     2 </A>
+            <A ID="31">   85     9    86    87 </A>
+            <A ID="32">   84    88    85    73 </A>
+            <A ID="33">   89    86    18    90 </A>
+            <A ID="34">   79    91    92    15 </A>
+            <A ID="35">   82    92    93    89 </A>
+            <A ID="36">   94    95    96    88 </A>
+            <A ID="37">   97    98    99   100 </A>
+            <A ID="38">   96   101    87    99 </A>
+            <A ID="39">  102   103   100   104 </A>
+            <A ID="40">  105   106   107    93 </A>
+            <A ID="41">  107   104    90   108 </A>
+            <A ID="42">  109   110    37   111 </A>
+            <A ID="43">  112    98   113   114 </A>
+            <A ID="44">  111   101   114    39 </A>
+            <A ID="45">  115   103   116   113 </A>
+            <A ID="46">  117   118    47   119 </A>
+            <A ID="47">  119   116   108    43 </A>
+        </ELEMENT>
+        <COMPOSITE>
+            <C ID="10"> A[0-10] </C>
+            <C ID="11"> A[11-25] </C>
+            <C ID="12"> A[26-40] </C>
+            <C ID="13"> A[41-47] </C>
+            <C ID="1"> F[0,19,34,49,64,83,94,109] </C>
+            <C ID="2"> F[1,4,10,13,65,68,74,78] </C>
+            <C ID="3"> F[69,72,75,81,95,97,102,105] </C>
+            <C ID="4"> F[50,52,57,61,110,112,115,117] </C>
+            <C ID="5"> F[20,23,27,29,53,56,58,63] </C>
+            <C ID="6"> F[14,30,45,60,77,91,106,118] </C>
+        </COMPOSITE>
+        <DOMAIN> C[10-13] </DOMAIN>
+    </GEOMETRY>
+    <EXPANSIONS>
+        <E COMPOSITE="C[10]" NUMMODES="4" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+        <E COMPOSITE="C[11]" NUMMODES="6" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+        <E COMPOSITE="C[12]" NUMMODES="4" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+        <E COMPOSITE="C[13]" NUMMODES="5" FIELDS="u,v,w,p" TYPE="MODIFIED" />
+    </EXPANSIONS>
+    <CONDITIONS>
+        <SOLVERINFO>
+            <I PROPERTY="SolverType" VALUE="VelocityCorrectionScheme" />
+            <I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes" />
+            <I PROPERTY="AdvectionForm" VALUE="Convective" />
+            <I PROPERTY="Projection" VALUE="Galerkin" />
+            <I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder1" />
+        </SOLVERINFO>
+
+        <PARAMETERS>
+            <P> TimeStep      = 0.001 </P>
+            <P> NumSteps      = 10 </P>
+            <P> IO_CheckSteps = 10    </P>
+            <P> IO_InfoSteps  = 10    </P>
+            <P> Kinvis        = 1   </P>
+        </PARAMETERS>
+
+        <VARIABLES>
+            <V ID="0"> u </V>
+            <V ID="1"> v </V>
+            <V ID="2"> w </V>
+            <V ID="3"> p </V>
+        </VARIABLES>
+
+        <BOUNDARYREGIONS>
+            <B ID="0"> C[1] </B>  <!-- Inlet -->
+            <B ID="1"> C[6] </B>    <!-- Outlet -->
+            <B ID="2"> C[2] </B>    <!-- Wall -->
+            <B ID="3"> C[3] </B>    <!-- Wall left -->
+            <B ID="4"> C[4] </B>    <!-- Wall -->
+            <B ID="5"> C[5] </B>    <!-- Wall right -->
+        </BOUNDARYREGIONS>
+
+        <BOUNDARYCONDITIONS>
+            <REGION REF="0">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="1">
+                <N VAR="u" VALUE="0" />
+                <N VAR="v" VALUE="0" />
+                <N VAR="w" VALUE="0" />
+                <D VAR="p" VALUE="0" />
+            </REGION>
+            <REGION REF="2">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="3">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="4">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="5">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="y*(1-y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+        </BOUNDARYCONDITIONS>
+
+        <FUNCTION NAME="ExactSolution">
+            <E VAR="u" VALUE="0" />
+            <E VAR="v" VALUE="0" />
+            <E VAR="w" VALUE="y*(1-y)" />
+            <E VAR="p" VALUE="-2*Kinvis*(z-1)" />
+        </FUNCTION>
+        <FUNCTION NAME="InitialConditions">
+            <E VAR="u" VALUE="0" />
+            <E VAR="v" VALUE="0" />
+            <E VAR="w" VALUE="y*(1-y)" />
+            <E VAR="p" VALUE="-2*Kinvis*(z-1)" />
+        </FUNCTION>
+    </CONDITIONS>
+</NEKTAR>
+