diff --git a/src/libvolumetricMesh/tetMesh.cpp b/src/libvolumetricMesh/tetMesh.cpp
index 8022ae8e4610d3d29a5ee576272b9773d37b9b8c..50615d82a69f175c0bb3de361bbe2996ee111c27 100644
--- a/src/libvolumetricMesh/tetMesh.cpp
+++ b/src/libvolumetricMesh/tetMesh.cpp
@@ -26,6 +26,8 @@
  *                                                                       *
  *************************************************************************/
 
+#include <cmath>
+
 #include "tetMesh.h"
 #include "volumetricMeshParser.h"
 
@@ -59,7 +61,7 @@ TetMesh::TetMesh(char * filename, int specialFileType, int verbose): VolumetricM
   sprintf(lineBuffer, "%s.node", filename);
   if (parser.open(lineBuffer) != 0)
     throw 2;
-  
+
   parser.getNextLine(lineBuffer, 1);
   int dim;
   sscanf(lineBuffer, "%d %d", &numVertices, &dim);
@@ -78,7 +80,7 @@ TetMesh::TetMesh(char * filename, int specialFileType, int verbose): VolumetricM
       throw 3;
     vertices[i] = new Vec3d(x,y,z);
   }
-  
+
   parser.close();
 
   // next, read the elements
@@ -109,10 +111,10 @@ TetMesh::TetMesh(char * filename, int specialFileType, int verbose): VolumetricM
     for(int j=0; j<4; j++) // vertices are 1-indexed in .ele files
     {
       v[j]--;
-      elements[i][j] = v[j]; 
+      elements[i][j] = v[j];
     }
   }
-  
+
   parser.close();
 
   numMaterials = 0;
@@ -133,7 +135,7 @@ TetMesh::TetMesh(int numVertices_, double * vertices_,
          int numElements_, int * elements_,
          int numMaterials_, Material ** materials_,
          int numSets_, Set ** sets_,
-         int numRegions_, Region ** regions_): 
+         int numRegions_, Region ** regions_):
   VolumetricMesh(numVertices_, vertices_, numElements_, 4, elements_,
                  numMaterials_, materials_, numSets_, sets_, numRegions_, regions_) {}
 
@@ -170,7 +172,7 @@ void TetMesh::computeElementMassMatrix(int el, double * massMatrix) const
    Singiresu S. Rao: The finite element method in engineering, 2004)
 */
 
-  const double mtx[16] = { 2, 1, 1, 1, 
+  const double mtx[16] = { 2, 1, 1, 1,
                            1, 2, 1, 1,
                            1, 1, 2, 1,
                            1, 1, 1, 2 } ;
@@ -308,7 +310,11 @@ double TetMesh::getTetDeterminant(Vec3d * a, Vec3d * b, Vec3d * c, Vec3d * d)
   Mat3d m2 = Mat3d(*a, *b, *d);
   Mat3d m3 = Mat3d(*a, *b, *c);
 
-  return (-det(m0) + det(m1) - det(m2) + det(m3));
+  double determinant = -det(m0) + det(m1) - det(m2) + det(m3);
+  if(std::fabs(determinant) < 1e-12)
+      determinant = 0.0;
+
+  return determinant;
 }
 
 void TetMesh::interpolateGradient(int element, const double * U, int numFields, Vec3d pos, double * grad) const
@@ -389,7 +395,7 @@ void TetMesh::getElementEdges(int el, int * edgeBuffer) const
     v[i] = getVertexIndex(el,i);
 
   int edgeMask[6][2] = {
-   { 0, 1 }, { 1, 2 }, { 2, 0 }, 
+   { 0, 1 }, { 1, 2 }, { 2, 0 },
    { 0, 3 }, { 1, 3 }, { 2, 3 } };
 
   for(int edge=0; edge<6; edge++)