Prismatic mesh

Code/Source/solver/VtkData.cpp

@@ -340,6 +340,9 @@ void VtkVtuData::VtkVtuDataImpl::set_connectivity(const int nsd, const Array<int
     } else if (np_elem == 4) {
       vtk_cell_type = VTK_TETRA;
 
+    } else if (np_elem == 6) {
+      vtk_cell_type = VTK_WEDGE;
+
     } else if (np_elem == 8) {
       vtk_cell_type = VTK_HEXAHEDRON;
 

Code/Source/solver/fs.cpp

@@ -302,13 +302,6 @@ void init_fs_msh(const ComMod& com_mod, mshType& lM)
     lM.fs[0].Nb  = lM.Nb;
     lM.fs[0].Nx  = lM.Nx;
   }
-
-  bool flag = (lM.eType == ElementType::HEX20) || (lM.eType == ElementType::HEX27) || (lM.eType == ElementType::WDG);
-  if (flag) {
-    std::cout << "WARNING: Second derivatives are not computed for HEX20/HEX27/WDG type elements";
-    //throw std::runtime_error(" Second derivatives are not computed for HEX20/HEX27/WDG type elements");
-  }
-
   // Second order derivatives for vector function space
   //
   if (!lM.fs[0].lShpF) {

Code/Source/solver/nn_elem_gip.h

@@ -265,6 +265,23 @@ GetElementGausIntMapType get_element_gauss_int_data = {
     xi(1,2) = s;
     }
   },
+
+  {ElementType::WDG, [](const int insd, const int nG, Vector<double>& w, Array<double>& xi)  -> void { 
+     for (int i = 0; i < nG; i++) { 
+       w(i) = 1.0 / 6.0;
+     }
+    double s  =  2.0 / 3.0;
+    double t  =  1.0 / 6.0;
+    double uz =  1.0 / sqrt(3.0);
+    double lz = -1.0 / sqrt(3.0);
+    xi(0,0) = s; xi(1,0) = t; xi(2,0) = lz;
+    xi(0,1) = t; xi(1,1) = s; xi(2,1) = lz;
+    xi(0,2) = t; xi(1,2) = t; xi(2,2) = lz;
+    xi(0,3) = s; xi(1,3) = t; xi(2,3) = uz;
+    xi(0,4) = t; xi(1,4) = s; xi(2,4) = uz;
+    xi(0,5) = t; xi(1,5) = t; xi(2,5) = uz;
+    }
+  },
 
 };
 

Code/Source/solver/nn_elem_props.h

@@ -61,7 +61,7 @@ SetElementPropsMapType set_3d_element_props = {
     mesh.nG = 6;
     mesh.vtkType = 13;
     mesh.nEf = 3;
-    mesh.lShpF = false;
+    mesh.lShpF = true;
     }
   },
 

Code/Source/solver/read_msh.cpp

@@ -758,12 +758,10 @@ void check_wedge_conn(mshType& mesh)
 
     if (sn == -1) {
       utils::swap( mesh.gIEN(0,e), mesh.gIEN(1,e));
-      a = 4;
-      b = 5;
-      if (e == 0) {
-        //std::cout << "[check_wedge_conn] Reorder element connectivity." << std::endl;
-      }
-    } else if (sn == 0) { 
+      a = 3;
+      b = 4;
+
+     } else if (sn == 0) { 
       throw std::runtime_error("Element " + std::to_string(e) + " is distorted.");
     }