diff --git a/Code/Source/svFSI/fluid.cpp b/Code/Source/svFSI/fluid.cpp
index e35b3362..2318e665 100644
--- a/Code/Source/svFSI/fluid.cpp
+++ b/Code/Source/svFSI/fluid.cpp
@@ -517,6 +517,9 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array<double>& Ag
   int num_c = lM.nEl / 10;
 
   for (int e = 0; e < lM.nEl; e++) {
+    #ifdef debug_construct_fluid
+    dmsg << "---------- e: " << e+1;
+    #endif
     cDmn = all_fun::domain(com_mod, lM, cEq, e);
     auto cPhys = eq.dmn[cDmn].phys;
 
@@ -579,7 +582,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array<double>& Ag
     //
     #ifdef debug_construct_fluid
     dmsg;
-    dmsg << "Gauss integration 1 ... ";
+    dmsg << "Gauss integration 1 ... " << "";
     dmsg << "fs[1].nG: " << fs[0].nG;
     dmsg << "fs[1].lShpF: " << fs[0].lShpF;
     dmsg << "fs[2].nG: " << fs[1].nG;
@@ -590,6 +593,9 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array<double>& Ag
     Array<double> ksix(nsd,nsd);
 
     for (int g = 0; g < fs[0].nG; g++) {
+      #ifdef debug_construct_fluid
+      dmsg << "===== g: " << g+1;
+      #endif
       if (g == 0 || !fs[1].lShpF) {
         auto Nx = fs[1].Nx.rslice(g);
         nn::gnn(fs[1].eNoN, nsd, nsd, Nx, xql, Nqx, Jac, ksix);
@@ -605,14 +611,15 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array<double>& Ag
            throw std::runtime_error("[construct_fluid] Jacobian for element " + std::to_string(e) + " is < 0.");
         }
 
-        if (!vmsStab) {
-          auto Nx = fs[0].Nx.rslice(g);
-          auto Nxx = fs[0].Nxx.rslice(g);
-          nn::gn_nxx(l, fs[0].eNoN, nsd, nsd, Nx, Nxx, xwl, Nwx, Nwxx); 
-        }
+        auto Nxx = fs[0].Nxx.rslice(g);
+        nn::gn_nxx(l, fs[0].eNoN, nsd, nsd, Nx, Nxx, xwl, Nwx, Nwxx); 
       }
 
       double w = fs[0].w(g) * Jac;
+      #ifdef debug_construct_fluid
+      dmsg << "Jac: " << Jac;
+      dmsg << "w: " << w;
+      #endif
 
       // Compute momentum residual and tangent matrix.
       //
@@ -636,6 +643,15 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array<double>& Ag
 
     // Gauss integration 2
     //
+    #ifdef debug_construct_fluid
+    dmsg;
+    dmsg << "Gauss integration 2 ... " << "";
+    dmsg << "fs[1].nG: " << fs[0].nG;
+    dmsg << "fs[1].lShpF: " << fs[0].lShpF;
+    dmsg << "fs[2].nG: " << fs[1].nG;
+    dmsg << "fs[2].lShpF: " << fs[1].lShpF;
+    #endif
+
     for (int g = 0; g < fs[1].nG; g++) {
       if (g == 0 || !fs[0].lShpF) {
         auto Nx = fs[0].Nx.rslice(g);
@@ -1205,7 +1221,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
   #endif
   
   // Maximum size of arrays sized by (3,eNoNw) -> (3,MAX_SIZE).
-  const int MAX_SIZE = 8;
+  const int MAX_SIZE = 27;
 
   using namespace consts;
 
@@ -1234,9 +1250,9 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
   // Velocity and its gradients, inertia (acceleration & body force)
   //
   double ud[3] = {-f[0], -f[1], -f[2]};
-  double u[3] = {0.0};
-  double ux[3][3] = {0.0};
-  double uxx[3][3][3] = {0.0};
+  double u[3] = {};
+  double ux[3][3] = {};
+  double uxx[3][3][3] = {};
 
   for (int a = 0; a < eNoNw; a++) {
     ud[0] = ud[0] + Nw(a)*(al(0,a)-bfl(0,a));
@@ -1257,7 +1273,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
     ux[1][2] = ux[1][2] + Nwx(1,a)*yl(2,a);
     ux[2][2] = ux[2][2] + Nwx(2,a)*yl(2,a);
 
-    uxx[0][0][1] += Nwxx(0,a)*yl(0,a);
+    uxx[0][0][0] += Nwxx(0,a)*yl(0,a);
     uxx[1][0][1] += Nwxx(1,a)*yl(0,a);
     uxx[2][0][2] += Nwxx(2,a)*yl(0,a);
     uxx[1][0][0] += Nwxx(3,a)*yl(0,a);
@@ -1293,14 +1309,14 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
   uxx[1][2][2] = uxx[2][2][1];
   uxx[2][2][0] = uxx[0][2][2];
 
-  double d2u2[3] = {0.0};
+  double d2u2[3] = {};
   d2u2[0] = uxx[0][0][0] + uxx[1][0][1] + uxx[2][0][2];
   d2u2[1] = uxx[0][1][0] + uxx[1][1][1] + uxx[2][1][2];
   d2u2[2] = uxx[0][2][0] + uxx[1][2][1] + uxx[2][2][2];
 
   // Pressure and its gradient
   //
-  double px[3] = {0.0};
+  double px[3] = {};
 
   for (int a = 0; a < eNoNq; a++) {
     px[0] = px[0] + Nqx(0,a)*yl(3,a);
@@ -1320,7 +1336,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
 
   // Strain rate tensor 2*e_ij := (u_i,j + u_j,i)
   //
-  double es[3][3] = {0.0};
+  double es[3][3] = {};
   es[0][0] = ux[0][0] + ux[0][0];
   es[1][1] = ux[1][1] + ux[1][1];
   es[2][2] = ux[2][2] + ux[2][2];
@@ -1339,7 +1355,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
     esNx[2][a] = es[0][2]*Nwx(0,a) + es[1][2]*Nwx(1,a) + es[2][2]*Nwx(2,a);
   }
 
-  double es_x[3][3][3] = {0.0};
+  double es_x[3][3][3] = {};
 
   for (int k = 0; k < 3; k++) { 
      es_x[0][0][k] = uxx[0][0][k] + uxx[0][0][k];
@@ -1391,13 +1407,11 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
   for (int i = 0; i < 3; i++) {
     mu_x[i] = mu_g * mu_x[i];
   }
-  //std::transform(mu_x.begin(), mu_x.end(), mu_x.begin(), [mu_g](double &v){return mu_g*v;});
-  //mu_x(:) = mu_g * mu_x(:)
 
   // Stabilization parameters
   //
-  double up[3] = {0.0};
-  double updu[3][3][MAX_SIZE] = {0.0};
+  double up[3] = {};
+  double updu[3][3][MAX_SIZE] = {};
   double tauM = 0.0;
 
   if (vmsFlag) {
@@ -1489,7 +1503,6 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
   }
 }
 
-
 /// @brief Element momentum residual.
 ///
 ///  Modifies:
@@ -1508,11 +1521,12 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
   dmsg << "vmsFlag: " << vmsFlag;
   dmsg << "eNoNw: " << eNoNw;
   dmsg << "eNoNq: " << eNoNq;
+  dmsg << "w: " << w;
   double start_time = utils::cput();
   #endif
 
   // Maximum size of arrays sized by (3,eNoNw) -> (3,MAX_SIZE).
-  const int MAX_SIZE = 8;
+  const int MAX_SIZE = 27;
 
   using namespace consts;
 
@@ -1549,9 +1563,9 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
   //
   std::array<double,3> ud{-f[0], -f[1], -f[2]};
   //ud  = -f
-  double u[3] = {0.0};
-  double ux[3][3] = {0.0};
-  double uxx[3][3][3] = {0.0};
+  double u[3] = {};
+  double ux[3][3] = {};
+  double uxx[3][3][3] = {};
 
   for (int a = 0; a < eNoNw; a++) {
     ud[0] = ud[0] + Nw(a)*(al(0,a)-bfl(0,a));
@@ -1619,7 +1633,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
   // Pressure and its gradient
   //
   double p = 0.0;
-  double px[3] = {0.0};
+  double px[3] = {};
 
   for (int a = 0; a < eNoNq; a++) {
     p  = p + Nq(a)*yl(3,a);
@@ -1628,6 +1642,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
     px[2] = px[2] + Nqx(2,a)*yl(3,a);
   }
   #ifdef debug_fluid_3d_m
+  dmsg << "u: " << u[0] << " " << u[1] << " " << u[2];
   dmsg << "p: " << p;
   dmsg << "px: " << px[0] << " " << px[1] << " " << px[2];
   #endif
@@ -1644,7 +1659,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
 
   // Strain rate tensor 2*e_ij := (u_i,j + u_j,i)
   //
-  double es[3][3] = {0.0};
+  double es[3][3] = {};
   es[0][0] = ux[0][0] + ux[0][0];
   es[1][1] = ux[1][1] + ux[1][1];
   es[2][2] = ux[2][2] + ux[2][2];
@@ -1745,24 +1760,24 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
   dmsg << "tauM: " << tauM;
   #endif
 
-  double rV[3] = {0.0};
+  double rV[3] = {};
   rV[0] = ud[0] + u[0]*ux[0][0] + u[1]*ux[1][0] + u[2]*ux[2][0];
   rV[1] = ud[1] + u[0]*ux[0][1] + u[1]*ux[1][1] + u[2]*ux[2][1];
   rV[2] = ud[2] + u[0]*ux[0][2] + u[1]*ux[1][2] + u[2]*ux[2][2];
 
-  double rS[3] = {0.0};
+  double rS[3] = {};
   rS[0] = mu_x[0]*es[0][0] + mu_x[1]*es[1][0] + mu_x[2]*es[2][0] + mu*d2u2[0];
   rS[1] = mu_x[0]*es[0][1] + mu_x[1]*es[1][1] + mu_x[2]*es[2][1] + mu*d2u2[1];
   rS[2] = mu_x[0]*es[0][2] + mu_x[1]*es[1][2] + mu_x[2]*es[2][2] + mu*d2u2[2];
 
-  double up[3] = {0.0};
+  double up[3] = {};
   up[0] = -tauM*(rho*rV[0] + px[0] - rS[0]);
   up[1] = -tauM*(rho*rV[1] + px[1] - rS[1]);
   up[2] = -tauM*(rho*rV[2] + px[2] - rS[2]);
 
   double tauC, tauB, pa;
   double eps = std::numeric_limits<double>::epsilon();
-  double ua[3] = {0.0};
+  double ua[3] = {};
 
   if (vmsFlag) {
     tauC = 1.0 / (tauM * (Kxi(0,0) + Kxi(1,1) + Kxi(2,2)));
@@ -1814,10 +1829,10 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
 
   //  Local residual
   //
-  double updu[3][3][MAX_SIZE] = {0.0};
-  double uNx[MAX_SIZE] = {0.0};
-  double upNx[MAX_SIZE] = {0.0}; 
-  double uaNx[MAX_SIZE] = {0.0}; 
+  double updu[3][3][MAX_SIZE] = {};
+  double uNx[MAX_SIZE] = {};
+  double upNx[MAX_SIZE] = {}; 
+  double uaNx[MAX_SIZE] = {}; 
 
   for (int a = 0; a < eNoNw; a++) {
     lR(0,a) = lR(0,a) + wr*Nw(a)*rV[0] + w*(Nwx(0,a)*rM[0][0] + Nwx(1,a)*rM[1][0] + Nwx(2,a)*rM[2][0]);
@@ -1825,6 +1840,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
     lR(2,a) = lR(2,a) + wr*Nw(a)*rV[2] + w*(Nwx(0,a)*rM[0][2] + Nwx(1,a)*rM[1][2] + Nwx(2,a)*rM[2][2]);
 
     // Quantities used for stiffness matrix
+
     uNx[a]  = u[0]*Nwx(0,a)  + u[1]*Nwx(1,a)  + u[2]*Nwx(2,a);
     upNx[a] = up[0]*Nwx(0,a) + up[1]*Nwx(1,a) + up[2]*Nwx(2,a);
 
@@ -1905,6 +1921,8 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
     }
   }
 
+  //exit(0);
+
   for (int b = 0; b < eNoNq; b++) {
     for (int a = 0; a < eNoNw; a++) {
       T1 = rho*tauM*uaNx[a];
diff --git a/Code/Source/svFSI/fsi.cpp b/Code/Source/svFSI/fsi.cpp
index b7010b93..1526f9cf 100644
--- a/Code/Source/svFSI/fsi.cpp
+++ b/Code/Source/svFSI/fsi.cpp
@@ -204,11 +204,8 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar
            throw std::runtime_error("[construct_fsi] Jacobian for element " + std::to_string(e) + " is < 0.");
         }
 
-        if (!vmsStab) {
-          auto Nx = fs_1[0].Nx.rslice(g);
-          auto Nxx = fs_1[0].Nxx.rslice(g);
-          nn::gn_nxx(l, fs_1[0].eNoN, nsd, nsd, Nx, Nxx, xwl, Nwx, Nwxx);
-        }
+        auto Nxx = fs_1[0].Nxx.rslice(g);
+        nn::gn_nxx(l, fs_1[0].eNoN, nsd, nsd, Nx, Nxx, xwl, Nwx, Nwxx);
       }
 
       double w = fs_1[0].w(g) * Jac;
diff --git a/Code/Source/svFSI/main.cpp b/Code/Source/svFSI/main.cpp
index c9eb5df4..d6ebd7ba 100644
--- a/Code/Source/svFSI/main.cpp
+++ b/Code/Source/svFSI/main.cpp
@@ -413,6 +413,7 @@ void iterate_solution(Simulation* simulation)
       #endif
 
       ls_ns::ls_alloc(com_mod, eq);
+      com_mod.Val.write("Val_alloc"+ istr);
 
       // Compute body forces. If phys is shells or CMM (init), apply
       // contribution from body forces (pressure) to residual
@@ -424,6 +425,7 @@ void iterate_solution(Simulation* simulation)
       #endif
 
       bf::set_bf(com_mod, Dg);
+      com_mod.Val.write("Val_bf"+ istr);
 
       // Assemble equations.
       //
@@ -437,9 +439,6 @@ void iterate_solution(Simulation* simulation)
       com_mod.R.write("R_as"+ istr);
       com_mod.Val.write("Val_as"+ istr);
 
-      com_mod.Val.write("Val_as"+ istr);
-      com_mod.R.write("R_as"+ istr);
-
       // Treatment of boundary conditions on faces
       // Apply Neumman or Traction boundary conditions
       //
diff --git a/Code/Source/svFSI/nn.cpp b/Code/Source/svFSI/nn.cpp
index 2c864f27..d0e5e862 100644
--- a/Code/Source/svFSI/nn.cpp
+++ b/Code/Source/svFSI/nn.cpp
@@ -815,7 +815,7 @@ void gn_nxx(const int l, const int eNoN, const int nsd, const int insd, Array<do
     }
 
     K.set_row(0, {xXi(0,0)*xXi(0,0), xXi(1,0)*xXi(1,0), t*xXi(0,0)*xXi(1,0)});
-    K.set_row(1, {xXi(0,1)*xXi(0,0), xXi(1,1)*xXi(1,1), t*xXi(0,1)*xXi(1,1)});
+    K.set_row(1, {xXi(0,1)*xXi(0,1), xXi(1,1)*xXi(1,1), t*xXi(0,1)*xXi(1,1)});
     K.set_row(2, {xXi(0,0)*xXi(0,1), xXi(1,0)*xXi(1,1), xXi(0,0)*xXi(1,1) + xXi(0,1)*xXi(1,0)});
 
     for (int a = 0; a < eNoN; a++) {
@@ -855,7 +855,7 @@ void gn_nxx(const int l, const int eNoN, const int nsd, const int insd, Array<do
     }
 
     for (int i = 0; i < 3; i++) { 
-      K.set_row(i, { xXi(0,i)*xXi(0,i), xXi(1,i)*xXi(1,i), xXi(2,i)*xXi(1,i), 
+      K.set_row(i, { xXi(0,i)*xXi(0,i), xXi(1,i)*xXi(1,i), xXi(2,i)*xXi(2,i), 
                      t*xXi(0,i)*xXi(1,i), t*xXi(1,i)*xXi(2,i), t*xXi(0,i)*xXi(2,i) 
                    } );
     }
@@ -885,7 +885,7 @@ void gn_nxx(const int l, const int eNoN, const int nsd, const int insd, Array<do
                     xXi(2,i)*xXi(2,j),
                     xXi(0,i)*xXi(1,j) + xXi(0,j)*xXi(1,i),
                     xXi(1,i)*xXi(2,j) + xXi(1,j)*xXi(2,i),
-                    xXi(0,i)*xXi(3,j) + xXi(0,j)*xXi(2,i) 
+                    xXi(0,i)*xXi(2,j) + xXi(0,j)*xXi(2,i) 
                   } );
 
 
diff --git a/Code/Source/svFSI/nn_elem_gnn.h b/Code/Source/svFSI/nn_elem_gnn.h
index d471911e..aa746ed2 100644
--- a/Code/Source/svFSI/nn_elem_gnn.h
+++ b/Code/Source/svFSI/nn_elem_gnn.h
@@ -1256,10 +1256,9 @@ SetElementShapeMapType set_element_shape_data = {
   {ElementType::TET10, [](int g, mshType& mesh) -> void {
     auto& xi = mesh.xi;
     auto& N = mesh.N;
-
     double s = 1.0 - xi(0,g) - xi(1,g) - xi(2,g);
     N(0,g)  = xi(0,g)*(2.0*xi(0,g) - 1.0);
-    N(1,g)  = xi(1,g)*(1.0*xi(1,g) - 1.0);
+    N(1,g)  = xi(1,g)*(2.0*xi(1,g) - 1.0);
     N(2,g)  = xi(2,g)*(2.0*xi(2,g) - 1.0);
     N(3,g)  = s    *(2.0*s     - 1.0);
     N(4,g)  = 4.0*xi(0,g)*xi(1,g);
diff --git a/Code/Source/svFSI/vtk_xml.cpp b/Code/Source/svFSI/vtk_xml.cpp
index 2d5804e0..27e87910 100644
--- a/Code/Source/svFSI/vtk_xml.cpp
+++ b/Code/Source/svFSI/vtk_xml.cpp
@@ -445,7 +445,16 @@ void read_vtp(const std::string& file_name, faceType& face)
     throw std::runtime_error("The VTK face file '" + file_name + "' can't be read.");
   }
 
-  vtk_xml_parser::load_vtp(file_name, face);
+  // Check if the face mesh needs to be read from a VTP or VTU file.
+  //
+  auto file_ext = file_name.substr(file_name.find_last_of(".") + 1);
+
+  if (file_ext == "vtp") {
+    vtk_xml_parser::load_vtp(file_name, face);
+
+  } else if (file_ext == "vtu") {
+    vtk_xml_parser::load_vtu(file_name, face);
+  }
 
   if (face.gN.size() == 0) {
     std::cout << "[WARNING] No node IDs found in the '" << file_name << "' face file.";
diff --git a/Code/Source/svFSI/vtk_xml_parser.cpp b/Code/Source/svFSI/vtk_xml_parser.cpp
index 3abd5eb6..74a77a1f 100644
--- a/Code/Source/svFSI/vtk_xml_parser.cpp
+++ b/Code/Source/svFSI/vtk_xml_parser.cpp
@@ -40,6 +40,7 @@
 #include <vtkDoubleArray.h>
 #include "vtkCellData.h"
 #include <vtkGenericCell.h>
+#include <vtkGeometryFilter.h>
 #include <vtkIntArray.h>
 #include <vtkPointData.h>
 #include <vtkDataArray.h>
@@ -143,104 +144,11 @@ const std::string ELEMENT_IDS_NAME("GlobalElementID");
 //             I n t e r n a l  U t i l i t i e s              //
 /////////////////////////////////////////////////////////////////
 
-/// @brief Store element connectivity into the Face object.
-///
-/// Face variables set
-///   face.nEl - number of elements
-///   face.eNoN - number of noders per element
-///   face.IEN - element connectivity (num_nodes_per_elem, num_elems)
+/// @brief Get the mesh nodes per element and ordering. 
 ///
-/// Note: The face.IEN array is allocated as in the Fortran code as (face.eNoN, face.nEl).
-//
-void store_element_conn(vtkSmartPointer<vtkPolyData> vtk_polydata, faceType& face)
-{
-  #define n_debug_store_element_conn
-  #ifdef debug_store_element_conn
-  std::cout << "[store_element_conn(polydata)] " << std::endl;
-  std::cout << "[store_element_conn(polydata)] ========== store_element_conn(polydata) =========" << std::endl;
-  #endif
-
-  vtkObject::GlobalWarningDisplayOff();
-
-  // Get the number of nodes per cell.
-  auto cell = vtkGenericCell::New();
-  vtk_polydata->GetCell(0, cell);
-  int np_elem = cell->GetNumberOfPoints();
-  #ifdef debug_store_element_conn
-  std::cout << "[store_element_conn(polydata)] np_elem: " << np_elem << std::endl;
-  std::cout << "[store_element_conn(polydata)] cell type: " << vtk_polydata->GetCellType(0) << std::endl;
-  #endif
-
-  // Allocate connectivity array.
-  auto num_elems = vtk_polydata->GetNumberOfCells();
-  face.nEl = num_elems;
-  face.eNoN = np_elem; 
-  face.IEN = Array<int>(np_elem, num_elems);
-
-  for (int i = 0; i < num_elems; i++) {
-    vtk_polydata->GetCell(i, cell);
-    auto num_cell_pts = cell->GetNumberOfPoints();
-    for (int j = 0; j < num_cell_pts; j++) {
-      auto id = cell->PointIds->GetId(j);
-      face.IEN(j,i) = id;
-    }
-  }
-}
-
-void store_element_conn(vtkSmartPointer<vtkPolyData> vtk_polydata, mshType& mesh)
+void get_mesh_ordering(const int num_elems, vtkSmartPointer<vtkUnsignedCharArray> cell_types, int& np_elem, 
+    std::vector<std::vector<int>>& ordering)
 {
-  #ifdef debug_store_element_conn
-  std::cout << "[store_element_conn(polydata,mesh)] " << std::endl;
-  std::cout << "[store_element_conn(polydata,mesh)] ========== store_element_conn(polydata,mesh) =========" << std::endl;
-  #endif
-
-  // Get the number of nodes per cell.
-  auto cell = vtkGenericCell::New();
-  vtk_polydata->GetCell(0, cell);
-  int np_elem = cell->GetNumberOfPoints();
-
-  // Allocate connectivity array.
-  auto num_elems = vtk_polydata->GetNumberOfCells();
-  mesh.gnEl = num_elems;
-  mesh.eNoN = np_elem;
-  mesh.gIEN = Array<int>(np_elem, num_elems);
-  #ifdef debug_store_element_conn
-  std::cout << "[store_element_conn(polydata,mesh)] num_elems: " << num_elems << std::endl;
-  std::cout << "[store_element_conn(polydata,mesh)] np_elem: " << np_elem << std::endl;
-  #endif
-
-  for (int i = 0; i < num_elems; i++) {
-    //int elem_id = elem_ids->GetValue(i);
-    vtk_polydata->GetCell(i, cell);
-    auto num_cell_pts = cell->GetNumberOfPoints();
-    for (int j = 0; j < num_cell_pts; j++) {
-      auto id = cell->PointIds->GetId(j);
-      mesh.gIEN(j,i) = id;
-    }
-  }
-}
-
-/// @brief Store VTK vtkUnstructuredGrid cell connectivity into a mesh element connectivity. 
-///
-/// Mesh variables set
-///   mesh.gnEl - number of elements
-///   mesh.eNoN - number of noders per element
-///   mesh.gIEN - element connectivity (num_nodes_per_elem, num_elems)
-///
-/// Note: The mesh.gIEN array is allocated as in the Fortran code as (np_elem, num_elems).
-//
-void store_element_conn(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, mshType& mesh)
-{
-  #ifdef debug_store_element_conn
-  std::cout << "[store_element_conn(ugrid)] " << std::endl;
-  std::cout << "[store_element_conn(ugrid)] ========== store_element_conn(ugrid) =========" << std::endl;
-  #endif
-  vtkSmartPointer<vtkUnsignedCharArray> cell_types = vtk_ugrid->GetCellTypesArray();
-  auto num_elems = vtk_ugrid->GetNumberOfCells();
-  #ifdef debug_store_element_conn
-  std::cout << "[store_element_conn(ugrid)] num_elems: " << num_elems << std::endl;
-  #endif
-
   int num_hex = 0;
   int num_line = 0;
   int num_quad = 0;
@@ -255,6 +163,7 @@ void store_element_conn(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, mshType&
   int num_quadratic_tetra = 0;
   int num_quadratic_hexahedron = 0;
   int num_triquadratic_hexahedron = 0;
+
   for (int i = 0; i < num_elems; i++) {
     switch (cell_types->GetValue(i)) {
       case VTK_HEXAHEDRON:
@@ -315,7 +224,7 @@ void store_element_conn(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, mshType&
     }
   }
 
-  #ifdef debug_store_element_conn
+  #ifdef debug_get_mesh_ordering
   std::cout << "[store_element_conn] num_line: " << num_line <<  std::endl;
   std::cout << "[store_element_conn] num_quad: " << num_quad <<  std::endl;
   std::cout << "[store_element_conn] num_hex: " << num_hex <<  std::endl;
@@ -324,8 +233,9 @@ void store_element_conn(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, mshType&
   std::cout << "[store_element_conn] num_unknown: " << num_unknown <<  std::endl;
   #endif
 
-  int np_elem = 0;
-  std::vector<std::vector<int>> ordering;
+  np_elem = 0;
+  ordering.clear();
+
   if (num_line != 0) {
     np_elem = vtk_cell_to_elem[VTK_LINE];
     ordering = vtk_cell_ordering[VTK_LINE];
@@ -370,6 +280,133 @@ void store_element_conn(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, mshType&
     np_elem = vtk_cell_to_elem[VTK_TRIQUADRATIC_HEXAHEDRON];
     ordering = vtk_cell_ordering[VTK_TRIQUADRATIC_HEXAHEDRON];
   }
+}
+
+/// @brief Store element connectivity into the Face object.
+///
+/// Face variables set
+///   face.nEl - number of elements
+///   face.eNoN - number of noders per element
+///   face.IEN - element connectivity (num_nodes_per_elem, num_elems)
+///
+/// Note: The face.IEN array is allocated as in the Fortran code as (face.eNoN, face.nEl).
+//
+void store_element_conn(vtkSmartPointer<vtkPolyData> vtk_polydata, faceType& face)
+{
+  #define n_debug_store_element_conn
+  #ifdef debug_store_element_conn
+  std::cout << "[store_element_conn(polydata)] " << std::endl;
+  std::cout << "[store_element_conn(polydata)] ========== store_element_conn(polydata) =========" << std::endl;
+  #endif
+
+  vtkObject::GlobalWarningDisplayOff();
+
+  // Get the number of nodes per cell.
+  auto cell = vtkGenericCell::New();
+  vtk_polydata->GetCell(0, cell);
+  int np_elem = cell->GetNumberOfPoints();
+  #ifdef debug_store_element_conn
+  std::cout << "[store_element_conn(polydata)] np_elem: " << np_elem << std::endl;
+  std::cout << "[store_element_conn(polydata)] cell type: " << vtk_polydata->GetCellType(0) << std::endl;
+  #endif
+
+  // Allocate connectivity array.
+  auto num_elems = vtk_polydata->GetNumberOfCells();
+  face.nEl = num_elems;
+  face.eNoN = np_elem; 
+  face.IEN = Array<int>(np_elem, num_elems);
+
+  for (int i = 0; i < num_elems; i++) {
+    vtk_polydata->GetCell(i, cell);
+    auto num_cell_pts = cell->GetNumberOfPoints();
+    for (int j = 0; j < num_cell_pts; j++) {
+      auto id = cell->PointIds->GetId(j);
+      face.IEN(j,i) = id;
+    }
+  }
+}
+
+void store_element_conn(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, faceType& face)
+{
+  // Get the number of nodes per cell.
+  auto cell = vtkGenericCell::New();
+  vtk_ugrid->GetCell(0, cell);
+  int np_elem = cell->GetNumberOfPoints();
+
+  auto num_elems = vtk_ugrid->GetNumberOfCells();
+  face.nEl = num_elems;
+  face.eNoN = np_elem; 
+  face.IEN = Array<int>(np_elem, num_elems);
+
+  for (int i = 0; i < num_elems; i++) {
+    vtk_ugrid->GetCell(i, cell);
+    auto num_cell_pts = cell->GetNumberOfPoints();
+    for (int j = 0; j < num_cell_pts; j++) {
+      auto id = cell->PointIds->GetId(j);
+      face.IEN(j,i) = id;
+    }
+  }
+}
+
+void store_element_conn(vtkSmartPointer<vtkPolyData> vtk_polydata, mshType& mesh)
+{
+  #ifdef debug_store_element_conn
+  std::cout << "[store_element_conn(polydata,mesh)] " << std::endl;
+  std::cout << "[store_element_conn(polydata,mesh)] ========== store_element_conn(polydata,mesh) =========" << std::endl;
+  #endif
+
+  // Get the number of nodes per cell.
+  auto cell = vtkGenericCell::New();
+  vtk_polydata->GetCell(0, cell);
+  int np_elem = cell->GetNumberOfPoints();
+
+  // Allocate connectivity array.
+  auto num_elems = vtk_polydata->GetNumberOfCells();
+  mesh.gnEl = num_elems;
+  mesh.eNoN = np_elem;
+  mesh.gIEN = Array<int>(np_elem, num_elems);
+  #ifdef debug_store_element_conn
+  std::cout << "[store_element_conn(polydata,mesh)] num_elems: " << num_elems << std::endl;
+  std::cout << "[store_element_conn(polydata,mesh)] np_elem: " << np_elem << std::endl;
+  #endif
+
+  for (int i = 0; i < num_elems; i++) {
+    //int elem_id = elem_ids->GetValue(i);
+    vtk_polydata->GetCell(i, cell);
+    auto num_cell_pts = cell->GetNumberOfPoints();
+    for (int j = 0; j < num_cell_pts; j++) {
+      auto id = cell->PointIds->GetId(j);
+      mesh.gIEN(j,i) = id;
+    }
+  }
+}
+
+/// @brief Store VTK vtkUnstructuredGrid cell connectivity into a mesh element connectivity. 
+///
+/// Mesh variables set
+///   mesh.gnEl - number of elements
+///   mesh.eNoN - number of noders per element
+///   mesh.gIEN - element connectivity (num_nodes_per_elem, num_elems)
+///
+/// Note: The mesh.gIEN array is allocated as in the Fortran code as (np_elem, num_elems).
+//
+void store_element_conn(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, mshType& mesh)
+{
+  #ifdef debug_store_element_conn
+  std::cout << "[store_element_conn(ugrid)] " << std::endl;
+  std::cout << "[store_element_conn(ugrid)] ========== store_element_conn(ugrid) =========" << std::endl;
+  #endif
+  vtkSmartPointer<vtkUnsignedCharArray> cell_types = vtk_ugrid->GetCellTypesArray();
+  auto num_elems = vtk_ugrid->GetNumberOfCells();
+  #ifdef debug_store_element_conn
+  std::cout << "[store_element_conn(ugrid)] num_elems: " << num_elems << std::endl;
+  #endif
+
+  // Get the mesh nodes per element and ordering. 
+  //
+  int np_elem = 0; 
+  std::vector<std::vector<int>> ordering;
+  get_mesh_ordering(num_elems, cell_types, np_elem, ordering);
 
   // For generic higher-order elements.
   //
@@ -445,6 +482,30 @@ void store_element_ids(vtkSmartPointer<vtkPolyData> vtk_polydata, faceType& face
   }
 }
 
+/// @brief Store element IDs into a faceType object.
+///
+/// Face data set
+///   face.gE
+//
+void store_element_ids(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, faceType& face)
+{
+  auto elem_ids = vtkIntArray::SafeDownCast(vtk_ugrid->GetCellData()->GetArray(ELEMENT_IDS_NAME.c_str()));
+  if (elem_ids == nullptr) {
+    throw std::runtime_error("No '" + ELEMENT_IDS_NAME + "' data of type Int32 found in VTK mesh.");
+    return;
+  }
+  #ifdef debug_store_element_ids
+  std::cout << "[store_element_ids] Allocate face.gE ... " <<  std::endl;
+  #endif
+  int num_elem_ids = elem_ids->GetNumberOfTuples();
+  face.gE = Vector<int>(num_elem_ids);
+  // [NOTE] It is not clear how these IDs are used but if they
+  // index into arrays or vectors then they need to be offset by -1.
+  for (int i = 0; i < num_elem_ids; i++) {
+    face.gE(i) = elem_ids->GetValue(i) - 1;
+  }
+}
+
 /// @brief Store nodal coordinates from the VTK points into a face.
 ///
 /// Face data set
@@ -507,6 +568,22 @@ void store_nodal_ids(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, mshType& me
   }
 }
 
+/// @brief Store VTK node IDs data array into a face nodal IDs.
+///
+void store_nodal_ids(vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid, faceType& face)
+{
+  vtkIdType num_nodes = vtk_ugrid->GetNumberOfPoints();
+  auto node_ids = vtkIntArray::SafeDownCast(vtk_ugrid->GetPointData()->GetArray(NODE_IDS_NAME.c_str()));
+  if (node_ids == nullptr) {
+    return;
+  }
+  face.gN = Vector<int>(num_nodes);
+  for (int i = 0; i < num_nodes; i++) {
+    face.gN(i) = node_ids->GetValue(i);
+    //std::cout << "[store_nodal_ids] mesh.gN(" << i << "): " << mesh.gN(i) << std::endl;
+  }
+}
+
 /// @brief Store VTK node IDs data array into a face nodal IDs.
 ///
 /// Note: This array appeats to be optional. 
@@ -621,6 +698,7 @@ void load_fiber_direction_vtu(const std::string& file_name, const std::string& d
 //
 void load_vtp(const std::string& file_name, faceType& face)
 {
+  #define n_debug_load_vtp 
   #ifdef debug_load_vtp 
   std::cout << "[load_vtp] " << std::endl;
   std::cout << "[load_vtp] ===== vtk_xml_parser.cpp::load_vtp ===== " << std::endl;
@@ -754,6 +832,52 @@ void load_vtu(const std::string& file_name, mshType& mesh)
   store_element_conn(vtk_ugrid, mesh);
 }
 
+/// @brief Store a surface mesh read from a VTK .vtu file into a Face object.
+//
+void load_vtu(const std::string& file_name, faceType& face)
+{
+  #define n_debug_load_vtu 
+  #ifdef debug_load_vtu 
+  std::cout << "[load_vtu] " << std::endl;
+  std::cout << "[load_vtu] ===== vtk_xml_parser::load_vtu ===== " << std::endl;
+  std::cout << "[load_vtu] file_name: " << file_name << std::endl;
+  #endif
+
+  auto reader = vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
+  reader->SetFileName(file_name.c_str());
+  reader->Update();
+  vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid = reader->GetOutput();
+
+  vtkIdType num_nodes = vtk_ugrid->GetNumberOfPoints();
+  if (num_nodes == 0) {
+    throw std::runtime_error("Failed reading the VTK file '" + file_name + "'.");
+  }
+
+  vtkIdType num_elems = vtk_ugrid->GetNumberOfCells();
+  auto cell = vtkGenericCell::New();
+  vtk_ugrid->GetCell(0, cell);
+  int np_elem = cell->GetNumberOfPoints();
+
+  #ifdef debug_load_vtu 
+  std::cout << "[load_vtu] Number of nodes: " << num_nodes << std::endl;
+  std::cout << "[load_vtu] Number of elements: " << num_elems << std::endl;
+  std::cout << "[load_vtu] Number of nodes per element: " << np_elem << std::endl;
+  #endif
+
+  // Store nodal coordinates.
+  auto points = vtk_ugrid->GetPoints();
+  store_nodal_coords(points, face);
+
+  // Store nodal IDs.
+  store_nodal_ids(vtk_ugrid, face);
+
+  // Store element connectivity.
+  store_element_conn(vtk_ugrid, face);
+
+  // Store element IDs.
+  store_element_ids(vtk_ugrid, face);
+}
+
 
 /// @brief Read a time series field from a VTK .vtu file.
 ///
diff --git a/Code/Source/svFSI/vtk_xml_parser.h b/Code/Source/svFSI/vtk_xml_parser.h
index dec6bf06..fa89d710 100644
--- a/Code/Source/svFSI/vtk_xml_parser.h
+++ b/Code/Source/svFSI/vtk_xml_parser.h
@@ -55,7 +55,10 @@ void load_vtp(const std::string& file_name, mshType& mesh);
 
 void load_vtu(const std::string& file_name, mshType& mesh);
 
+void load_vtu(const std::string& file_name, faceType& face);
+
 void load_time_varying_field_vtu(const std::string file_name, const std::string field_name, mshType& mesh);
+
 };
 
 #endif
diff --git a/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-complete.mesh.vtu b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-complete.mesh.vtu
new file mode 100644
index 00000000..43f38266
--- /dev/null
+++ b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-complete.mesh.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:612fbff8ccfbe3e4616987c235a89075fe1e57653193c02271917caea1769e95
+size 38993
diff --git a/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/RPA6.vtu b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/RPA6.vtu
new file mode 100644
index 00000000..2626db0c
--- /dev/null
+++ b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/RPA6.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6cf7879a01648c677ef183a12db48fe69ae47970b20270b1bfd2025bc0684351
+size 1945
diff --git a/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/RPA61.vtu b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/RPA61.vtu
new file mode 100644
index 00000000..5d005133
--- /dev/null
+++ b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/RPA61.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:260aaeae4c42d41908d0bc22f8dd6c4e031e6160c4d817ba1428c3934c8779af
+size 1949
diff --git a/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/inflow.vtu b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/inflow.vtu
new file mode 100644
index 00000000..2998c43e
--- /dev/null
+++ b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/inflow.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5e1c39d1ab1d3e8f2cc85e3467130ca7a757cc19e9bcd5613ef2f98826b725da
+size 2450
diff --git a/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/wall_RPA6.vtu b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/wall_RPA6.vtu
new file mode 100644
index 00000000..14a5ea7f
--- /dev/null
+++ b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/wall_RPA6.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:40e5790b6162fa5bc21c074f7ba0336a75a3337187ce4393d9beba327779f0c3
+size 9037
diff --git a/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/wall_RPA61.vtu b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/wall_RPA61.vtu
new file mode 100644
index 00000000..b6a1ffe1
--- /dev/null
+++ b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/wall_RPA61.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a02016bde1e5ddde5ffdd178a03ae1a4b49ceb990322268ef827d5d333ccd612
+size 4987
diff --git a/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/wall_blend_RPA6_RPA61.vtu b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/wall_blend_RPA6_RPA61.vtu
new file mode 100644
index 00000000..69610af4
--- /dev/null
+++ b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/mesh-surfaces/wall_blend_RPA6_RPA61.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fd8eb036fa54b8918712966f2877ce0953dbb24cd14b8e4acd3f945ceff6e785
+size 2361
diff --git a/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/walls_combined.vtu b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/walls_combined.vtu
new file mode 100644
index 00000000..5ddb8022
--- /dev/null
+++ b/tests/cases/fluid/quadratic_tet10/quad-mesh-complete/walls_combined.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b241700259664ea79c72ceb77b9f7721cbc6fae18306733fb9a736bd178ab637
+size 12607
diff --git a/tests/cases/fluid/quadratic_tet10/result_003.vtu b/tests/cases/fluid/quadratic_tet10/result_003.vtu
new file mode 100644
index 00000000..4261afc3
--- /dev/null
+++ b/tests/cases/fluid/quadratic_tet10/result_003.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:21469b25ceeeb7ecead0d66458ca42766a1f4c385c534f4445afbb71d04fe46f
+size 224689
diff --git a/tests/cases/fluid/quadratic_tet10/svFSI.xml b/tests/cases/fluid/quadratic_tet10/svFSI.xml
new file mode 100755
index 00000000..0be4a7b6
--- /dev/null
+++ b/tests/cases/fluid/quadratic_tet10/svFSI.xml
@@ -0,0 +1,106 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<svFSIFile version="0.1">
+<GeneralSimulationParameters>
+
+  <Continue_previous_simulation> 0 </Continue_previous_simulation>
+  <Number_of_spatial_dimensions> 3 </Number_of_spatial_dimensions>
+  <Number_of_time_steps> 3 </Number_of_time_steps>
+  <Time_step_size> 0.005 </Time_step_size>
+  <Spectral_radius_of_infinite_time_step> 0.5 </Spectral_radius_of_infinite_time_step>
+  <Searched_file_name_to_trigger_stop> STOP_SIM </Searched_file_name_to_trigger_stop>
+
+  <Save_results_to_VTK_format> 1 </Save_results_to_VTK_format>
+  <Name_prefix_of_saved_VTK_files> result </Name_prefix_of_saved_VTK_files>
+  <Increment_in_saving_VTK_files> 3 </Increment_in_saving_VTK_files>
+  <Start_saving_after_time_step> 1 </Start_saving_after_time_step>
+
+  <Increment_in_saving_restart_files> 3 </Increment_in_saving_restart_files>
+  <Convert_BIN_to_VTK_format> 0 </Convert_BIN_to_VTK_format>
+
+  <Verbose> 1 </Verbose>
+  <Warning> 0 </Warning>
+  <Debug> 0 </Debug>
+
+</GeneralSimulationParameters>
+
+<Add_mesh name="msh" >
+    <Mesh_file_path> quad-mesh-complete/mesh-complete.mesh.vtu </Mesh_file_path>
+    <Add_face name="inflow">
+        <Face_file_path> quad-mesh-complete/mesh-surfaces/inflow.vtu </Face_file_path>
+    </Add_face>
+    <Add_face name="outlet1">
+        <Face_file_path> quad-mesh-complete/mesh-surfaces/RPA6.vtu </Face_file_path>
+    </Add_face>
+    <Add_face name="outlet2">
+        <Face_file_path> quad-mesh-complete/mesh-surfaces/RPA61.vtu </Face_file_path>
+    </Add_face>
+    <Add_face name="wall">
+        <Face_file_path> quad-mesh-complete/walls_combined.vtu </Face_file_path>
+    </Add_face>
+</Add_mesh>
+
+
+<Add_equation type="fluid" >
+   <Coupled> true </Coupled>
+   <Min_iterations> 3 </Min_iterations> 
+   <Max_iterations> 12 </Max_iterations> 
+   <Tolerance> 1e-7 </Tolerance> 
+   <Backflow_stabilization_coefficient> 0.2 </Backflow_stabilization_coefficient>
+   <Density> 1.06 </Density>
+
+   <Viscosity model="Constant" >
+     <Value> 0.04 </Value>
+   </Viscosity>
+
+   <Output type="Spatial" >
+     <Velocity> true </Velocity>
+     <Pressure> true </Pressure>
+     <Traction> true </Traction>
+     <WSS> true </WSS>
+     <Vorticity> true </Vorticity>
+     <Divergence> true </Divergence>
+
+   </Output>
+
+   <LS type="NS" >
+      <Linear_algebra type="fsils" >
+         <Preconditioner> fsils </Preconditioner>
+      </Linear_algebra> 
+      <Max_iterations> 15 </Max_iterations>
+      <NS_GM_max_iterations> 10 </NS_GM_max_iterations>
+      <NS_CG_max_iterations> 300 </NS_CG_max_iterations>
+      <Tolerance> 1e-3 </Tolerance>
+      <NS_GM_tolerance> 1e-3 </NS_GM_tolerance>
+      <NS_CG_tolerance> 1e-3 </NS_CG_tolerance>
+   </LS>
+
+   <Add_BC name="inflow" > 
+      <Type> Dir </Type>
+      <Time_dependence> Steady </Time_dependence>
+      <Value> -10.0 </Value>
+      <Impose_flux> true </Impose_flux>
+   </Add_BC>
+
+   <Add_BC name="outlet1" >
+      <Type> Neu </Type>
+      <Time_dependence> Resistance </Time_dependence>
+      <Value> 2666.0 </Value>
+   </Add_BC>
+
+
+  <Add_BC name="outlet2" >
+      <Type> Neu </Type>
+      <Time_dependence> Resistance </Time_dependence>
+      <Value> 2666.0 </Value>
+   </Add_BC>
+
+
+   <Add_BC name="wall" >
+      <Type> Dir </Type>
+      <Time_dependence> Steady </Time_dependence>
+      <Value> 0.0 </Value>
+   </Add_BC>
+
+</Add_equation>
+
+</svFSIFile>
diff --git a/tests/test_fluid.py b/tests/test_fluid.py
index 579bb7e3..39ed503c 100644
--- a/tests/test_fluid.py
+++ b/tests/test_fluid.py
@@ -86,3 +86,10 @@ def test_carreau_yasuda(n_proc):
 def test_iliac_artery(n_proc):
     test_folder = "iliac_artery"
     run_with_reference(base_folder, test_folder, fields, n_proc)
+
+
+def test_quadratic_tet10(n_proc):
+    test_folder = "quadratic_tet10"
+    t_max = 3
+    fields = ["Velocity", "Pressure", "Vorticity", "Divergence"]
+    run_with_reference(base_folder, test_folder, fields, n_proc, t_max)