diff --git a/libsrc/meshing/debugging.cpp b/libsrc/meshing/debugging.cpp
index feba17ab..c29aa688 100644
--- a/libsrc/meshing/debugging.cpp
+++ b/libsrc/meshing/debugging.cpp
@@ -97,7 +97,7 @@ namespace netgen
         return mesh_ptr;
     }
 
-    void CheckMesh (const Mesh& mesh, MESHING_STEP step)
+    void CheckMesh (const Mesh& mesh, MESHING_STEP step, const char * filename, int line)
     {
       if (step == MESHCONST_OPTVOLUME)
         {
@@ -107,12 +107,24 @@ namespace netgen
               double volume = el.Volume(mesh.Points());
               if (volume < 0)
                 {
+                  if(!have_error && line != -1)
+                    cerr << "Negative volume in mesh at " << filename << ":" << line << endl;
                   have_error = true;
-                  cout << "volume of element " << el << " is negative: " << volume << endl;
+                  cerr << "volume of element " << el << " is negative: " << volume << endl;
                 }
             }
           if (have_error)
-            throw Exception("Negative volume");
+            {
+              string s;
+              if(line != -1)
+                {
+                  s += filename;
+                  s += ":";
+                  s += ToString(line);
+                  s += "\t";
+                }
+              throw Exception(s + "Negative volume");
+            }
 
         CheckElementsAroundEdges(mesh);
         }
diff --git a/libsrc/meshing/debugging.hpp b/libsrc/meshing/debugging.hpp
index 722fffa0..d67f0e04 100644
--- a/libsrc/meshing/debugging.hpp
+++ b/libsrc/meshing/debugging.hpp
@@ -8,7 +8,7 @@ namespace netgen
     unique_ptr<Mesh> FilterMesh( const Mesh & m, FlatArray<PointIndex> points, FlatArray<SurfaceElementIndex> sels = Array<SurfaceElementIndex>{}, FlatArray<ElementIndex> els = Array<ElementIndex>{} );
 
     // Checks if the mesh is valid. This is called automatically on various places if debugparam.slowchecks is set
-    void CheckMesh( const Mesh & m, MESHING_STEP meshing_step );
+    void CheckMesh( const Mesh & m, MESHING_STEP meshing_step, const char *filename = "", int line = -1 );
 
     // Sometimes during SwapImprove we discover topological errors in the mesh. For instance, an edge is adjacent to 8 tets around it, but
     // the 8 "other" points of the tets don't form a closed path around the edge. Instead there are 2 sets of 4 points/tets each, which are not connected.
diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp
index d81e7a81..65f65e65 100644
--- a/libsrc/meshing/meshfunc.cpp
+++ b/libsrc/meshing/meshfunc.cpp
@@ -704,6 +704,7 @@ namespace netgen
     */
 
     mesh3d.CalcSurfacesOfNode();
+    CheckMesh(mesh3d, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
 
     MeshOptimize3d optmesh(mesh3d, mp);
 
@@ -714,6 +715,7 @@ namespace netgen
     bool do_swap2 = mp.optimize3d.find('t') != string::npos;
     for([[maybe_unused]] auto i : Range(mp.optsteps3d))
       {
+        CheckMesh(mesh3d, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
         auto [total_badness, max_badness, bad_els] = optmesh.UpdateBadness();
         if(bad_els==0) break;
         if(do_split) optmesh.SplitImprove();
@@ -724,6 +726,8 @@ namespace netgen
     // Now optimize all elements
     optmesh.SetMinBadness(0);
 
+    CheckMesh(mesh3d, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
+
     for (auto i : Range(mp.optsteps3d))
       {
 	if (multithread.terminate)
@@ -737,6 +741,8 @@ namespace netgen
 	    if (multithread.terminate)
 	      break;
 
+      CheckMesh(mesh3d, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
+
 	    switch (mp.optimize3d[j])
 	      {
 	      case 'c': 
@@ -762,6 +768,8 @@ namespace netgen
 	// mesh3d.mglevels = 1;
 	MeshQuality3d (mesh3d);
       }
+
+    CheckMesh(mesh3d, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
   
     multithread.task = savetask;
     return MESHING3_OK;
@@ -877,12 +885,19 @@ namespace netgen
       MeshOptimize3d optmesh(mesh, dummymp, OPT_CONFORM);
 
       for ([[maybe_unused]] auto i : Range(3)) {
+        CheckMesh(mesh, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
         optmesh.ImproveMesh();
+        CheckMesh(mesh, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
         optmesh.SwapImprove2 ();
+        CheckMesh(mesh, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
         optmesh.ImproveMesh();
+        CheckMesh(mesh, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
         optmesh.SwapImprove();
+        CheckMesh(mesh, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
         optmesh.ImproveMesh();
+        CheckMesh(mesh, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
         optmesh.CombineImprove();
+        CheckMesh(mesh, MESHCONST_OPTVOLUME, __FILE__, __LINE__);
       }
       last_num_bad_segs = num_bad_segs;
     }