High order element face mesh fails 205

Code/Source/svFSI/fluid.cpp

@@ -517,6 +517,9 @@
   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 @@
     //
     #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 @@
     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 @@
            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 @@
 
     // 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);
@@ -1217,7 +1233,7 @@
   #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;
 
@@ -1246,9 +1262,9 @@
   // 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));
@@ -1269,7 +1285,7 @@
     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);
@@ -1305,14 +1321,14 @@
   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);
@@ -1332,7 +1348,7 @@
 
   // 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];
@@ -1351,7 +1367,7 @@
     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];
@@ -1403,13 +1419,11 @@
   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) {
@@ -1501,7 +1515,6 @@
   }
 }
 
-
 /// @brief Element momentum residual.
 ///
 ///  Modifies:
@@ -1520,11 +1533,12 @@
   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;
 
@@ -1561,9 +1575,9 @@
   //
   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));
@@ -1631,7 +1645,7 @@
   // 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);
@@ -1640,6 +1654,7 @@
     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
@@ -1656,7 +1671,7 @@
 
   // 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];
@@ -1757,24 +1772,24 @@
   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)));
@@ -1826,17 +1841,18 @@
 
   //  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]);
     lR(1,a) = lR(1,a) + wr*Nw(a)*rV[1] + w*(Nwx(0,a)*rM[0][1] + Nwx(1,a)*rM[1][1] + Nwx(2,a)*rM[2][1]);
     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);
 
@@ -1917,6 +1933,8 @@
     }
   }
 
+  //exit(0);
+
   for (int b = 0; b < eNoNq; b++) {
     for (int a = 0; a < eNoNw; a++) {
       T1 = rho*tauM*uaNx[a];

Code/Source/svFSI/fsi.cpp

@@ -207,11 +207,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;

Code/Source/svFSI/main.cpp

@@ -308,6 +308,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
@@ -319,6 +320,7 @@ void iterate_solution(Simulation* simulation)
       #endif
 
       bf::set_bf(com_mod, Dg);
+      com_mod.Val.write("Val_bf"+ istr);
 
       // Assemble equations.
       //
@@ -332,9 +334,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
       //

Code/Source/svFSI/nn.cpp

@@ -788,7 +788,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++) {
@@ -828,7 +828,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) 
                    } );
     }
@@ -858,7 +858,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) 
                   } );
 
 

Code/Source/svFSI/nn_elem_gnn.h

@@ -1256,10 +1256,9 @@
   {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);

Code/Source/svFSI/vtk_xml.cpp

@@ -445,7 +445,16 @@
     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.";