diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index f2d0e0f96fb857fb5e95cfaeedade227b0fa465c..a20f880db31b8bbc3f4a900fb255a03ac6fcf1a5 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -671,6 +671,9 @@ static void Mesh3D(GModel *m)
     }
   }
 
+  // Ensure that all volume Jacobians are positive
+  m->setAllVolumesPositive();
+
   double t2 = Cpu();
   CTX::instance()->meshTimer[2] = t2 - t1;
   Msg::StatusBar(true, "Done meshing 3D (%g s)", CTX::instance()->meshTimer[2]);
@@ -683,6 +686,9 @@ void OptimizeMeshNetgen(GModel *m)
 
   std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegionNetgen());
 
+  // Ensure that all volume Jacobians are positive
+  m->setAllVolumesPositive();
+
   double t2 = Cpu();
   Msg::StatusBar(true, "Done optimizing 3D mesh with Netgen (%g s)", t2 - t1);
 }
@@ -694,6 +700,9 @@ void OptimizeMesh(GModel *m)
 
   std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegionGmsh());
 
+  // Ensure that all volume Jacobians are positive
+  m->setAllVolumesPositive();
+
   double t2 = Cpu();
   Msg::StatusBar(true, "Done optimizing 3D mesh (%g s)", t2 - t1);
 }
@@ -782,8 +791,6 @@ void GenerateMesh(GModel *m, int ask)
     }
   }
 
-  m->setAllVolumesPositive();
-
   // Subdivide into quads or hexas
   if(m->getMeshStatus() == 2 && CTX::instance()->mesh.algoSubdivide == 1)
     RefineMesh(m, CTX::instance()->mesh.secondOrderLinear, true);