Fiber direction spatial output seg fault 432

Code/Source/solver/post.cpp

@@ -681,27 +681,31 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra
       int Ac = lM.IEN(a,e);
       for (int i = 0; i < nsd; i++) {
         xl(i,a) = com_mod.x(i,Ac);
+      }
+      for (int i = 0; i < tDof; i++) {
         dl(i,a) = lD(i,Ac);
       }
     }
 
     for (int iFn = 0; iFn < lM.nFn; iFn++) {
       for (int i = 0; i < nsd; i++) {
-        fN(i,iFn) = lM.fN(iFn*nsd,e);
+        fN(i,iFn) = lM.fN(i+iFn*nsd,e);
       }
     }
 
-    Array<double> F;
+    Array<double> F(nsd,nsd);
+    Array<double> Nx(nsd,eNoN);
+    double Jac = 0.0;
 
     for (int g = 0; g < lM.nG; g++) {
-      double Jac = 0.0;
       if (g == 0 || !lM.lShpF) {
-        auto Nx = lM.Nx.slice(g);
-        nn::gnn(eNoN, nsd, nsd, Nx, xl, Nx, Jac, F);
+        auto Nxi = lM.Nx.slice(g);
+        nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, F);
       }
 
       double w = lM.w(g) * Jac;
-      auto F = mat_fun::mat_id(nsd); 
+      N = lM.N.col(g);
+      F = mat_fun::mat_id(nsd); 
 
       for (int a = 0; a < eNoN; a++) {
         if (nsd == 3) {
@@ -749,11 +753,12 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra
   for (int a = 0; a < lM.nNo; a++) {
     int Ac = lM.gN(a);
     if (!utils::is_zero(sA(Ac))) {
-      for (int i = 0; i < nsd; i++) {
+      for (int i = 0; i < nFn*nsd; i++) {
         res(i,a) = res(i,a) + sF(i,Ac) / sA(Ac);
-       }
+      }
     } 
   }
+
 }
 
 /// @brief Compute fiber stretch based on 4th invariant: I_{4,f}

Code/Source/solver/vtk_xml.cpp

@@ -86,8 +86,6 @@ void do_test()
     }
   }
 
-  std::cout << "Num dupe: " << num_dupe << std::endl;
-
   vtkSmartPointer<vtkUnstructuredGrid> mesh;
 
   std::string fileName = "mesh-complete/mesh-complete.mesh.vtu";
@@ -98,8 +96,6 @@ void do_test()
 
   vtkIdType m_NumPoints = mesh->GetNumberOfPoints();
   vtkIdType m_NumCells = mesh->GetNumberOfCells();
-  std::cout << "  Number of points: " << m_NumPoints << std::endl;
-  std::cout << "  Number of cells: " << m_NumCells << std::endl;
 
   auto numPoints = mesh->GetNumberOfPoints();
   auto points = mesh->GetPoints();
@@ -138,10 +134,6 @@ void do_test()
     }
   }
 
-  std::cout << "num_found_1: " << num_found_1 << std::endl;
-  std::cout << "num_found_2: " << num_found_2 << std::endl;
-
-  //exit(0);
 }
 
 /// @brief This routine prepares data array of a regular mesh
@@ -1171,11 +1163,11 @@ void write_vtus(Simulation* simulation, const Array<double>& lA, const Array<dou
               post::fib_dir_post(simulation, msh, nFn, tmpV, lD, iEq);
             }
             for (int iFn = 0; iFn < nFn; iFn++) {
+              cOut = cOut + 1;
               is = outS[cOut];
               ie = is + l - 1;
               outS[cOut+1] = ie + 1;
-              outNames[cOut] = eq.output[iOut].name + std::to_string(iFn);
-              cOut = cOut + 1;
+              outNames[cOut] = eq.output[iOut].name + std::to_string(iFn+1);
 
               for (int a = 0; a < msh.nNo; a++) {
                 for (int i = 0; i < l; i++) {
@@ -1414,8 +1406,6 @@ void write_vtus(Simulation* simulation, const Array<double>& lA, const Array<dou
     nSh = nSh + d[iM].nNo;
   }
 
-  // std::cout << "nOut: " << nOut << std::endl;
-
   // Writing all solutions
   //
   for (int iOut = 1; iOut < nOut; iOut++) {
@@ -1436,9 +1426,6 @@ void write_vtus(Simulation* simulation, const Array<double>& lA, const Array<dou
       nSh = nSh + d[iM].nNo;
     }
 
-    // std::cout << "iOut: " << iOut << std::endl;
-    // std::cout << "outName: " << outNames[iOut] << std::endl;
-
     vtk_writer->set_point_data(outNames[iOut], tmpV);
   }