diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index c75fa1f80d50e11e6ddab6dc9b60df2ece62859a..5a888190841e7886d4548a5405dc3126223318ac 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -1,40 +1,29 @@
 #ifndef CONNECTIVITY_HPP
 #define CONNECTIVITY_HPP
 
-#include <utils/PugsAssert.hpp>
-#include <utils/PugsMacros.hpp>
-
-#include <utils/PugsUtils.hpp>
-
-#include <utils/PugsTraits.hpp>
-
 #include <algebra/TinyVector.hpp>
-
-#include <mesh/ItemValue.hpp>
-
-#include <mesh/IConnectivity.hpp>
-
+#include <mesh/CellType.hpp>
+#include <mesh/ConnectivityComputer.hpp>
 #include <mesh/ConnectivityMatrix.hpp>
+#include <mesh/IConnectivity.hpp>
 #include <mesh/ItemToItemMatrix.hpp>
-
-#include <mesh/ConnectivityComputer.hpp>
-#include <mesh/SubItemValuePerItem.hpp>
-
-#include <algorithm>
-#include <vector>
-
-#include <mesh/CellType.hpp>
-
-#include <utils/CSRGraph.hpp>
-
 #include <mesh/ItemType.hpp>
+#include <mesh/ItemValue.hpp>
 #include <mesh/RefId.hpp>
 #include <mesh/RefItemList.hpp>
-
+#include <mesh/SubItemValuePerItem.hpp>
 #include <mesh/SynchronizerManager.hpp>
+#include <utils/CSRGraph.hpp>
+#include <utils/Exceptions.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+#include <utils/PugsTraits.hpp>
+#include <utils/PugsUtils.hpp>
 
+#include <algorithm>
 #include <iostream>
 #include <set>
+#include <vector>
 
 class ConnectivityDescriptor;
 
@@ -208,8 +197,7 @@ class Connectivity final : public IConnectivity
   const auto&
   edgeOwner() const
   {
-    std::cerr << __FILE__ << ':' << __LINE__ << ": edge owner not built\n";
-    std::terminate();
+    throw NotImplementedError("edge owner not built");
     return m_edge_owner;
   }
 
@@ -255,8 +243,7 @@ class Connectivity final : public IConnectivity
   const auto&
   edgeIsOwned() const
   {
-    std::cerr << __FILE__ << ':' << __LINE__ << ": edge is owned not built\n";
-    std::terminate();
+    throw NotImplementedError("edge is owned not built");
     return m_edge_is_owned;
   }
 
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index f74316c566dc69992eddc7bf5dfabdf537be5195..46deb5b95166d42bc3bb5ade5d3a02683157abfa 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -1,10 +1,11 @@
-#include <mesh/Connectivity.hpp>
-
 #include <mesh/ConnectivityComputer.hpp>
+
+#include <mesh/Connectivity.hpp>
 #include <utils/Exceptions.hpp>
 
 #include <iostream>
 #include <map>
+#include <sstream>
 
 template <typename ConnectivityType>
 PUGS_INLINE ConnectivityMatrix
@@ -20,9 +21,9 @@ ConnectivityComputer::computeConnectivityMatrix(const ConnectivityType& connecti
 
     item_to_child_item_matrix = this->_computeInverse(child_to_item_matrix);
   } else {
-    std::cerr << "unable to compute connectivity " << itemName(item_type) << " -> " << itemName(child_item_type)
-              << '\n';
-    std::terminate();
+    std::stringstream error_msg;
+    error_msg << "unable to compute connectivity " << itemName(item_type) << " -> " << itemName(child_item_type);
+    throw UnexpectedError(error_msg.str());
   }
 
   return item_to_child_item_matrix;
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 0d38e108bda1fdfa2be5bc5d3a742648ca3ccdc6..629f32801454edcf3e25a3c68aeeaa972a4ea6ff 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -643,8 +643,9 @@ GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& d
         break;
       }
       default: {
-        std::cerr << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 2!\n";
-        std::terminate();
+        std::ostringstream error_msg;
+        error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 2";
+        throw UnexpectedError(error_msg.str());
       }
       }
     } else if constexpr (Dimension == 3) {
@@ -697,8 +698,9 @@ GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& d
         break;
       }
       default: {
-        std::cerr << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3!\n";
-        std::terminate();
+        std::ostringstream error_msg;
+        error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
+        throw UnexpectedError(error_msg.str());
       }
       }
     }
@@ -843,8 +845,7 @@ GmshReader::__computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDe
           Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
           auto i = edge_id_map.find(edge);
           if (i == edge_id_map.end()) {
-            std::cerr << "could not find this edge!\n";
-            std::terminate();
+            throw NormalError("could not find this edge");
           }
           cell_edge_vector.push_back(i->second);
         }
@@ -861,8 +862,7 @@ GmshReader::__computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDe
           Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
           auto i = edge_id_map.find(edge);
           if (i == edge_id_map.end()) {
-            std::cerr << "could not find this edge!\n";
-            std::terminate();
+            throw NormalError("could not find this edge");
           }
           cell_edge_vector.push_back(i->second);
         }
@@ -870,8 +870,9 @@ GmshReader::__computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDe
         break;
       }
       default: {
-        std::cerr << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3!\n";
-        std::terminate();
+        std::stringstream error_msg;
+        error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
+        throw UnexpectedError(error_msg.str());
       }
       }
     }
@@ -1206,8 +1207,7 @@ GmshReader::__proceedData()
           auto i =
             face_to_id_map.find(Face({__triangles[f][0], __triangles[f][1], __triangles[f][2]}, node_number_vector));
           if (i == face_to_id_map.end()) {
-            std::cerr << "face not found!\n";
-            std::terminate();
+            throw NormalError("face not found");
           }
           return i->second;
         }();
@@ -1239,8 +1239,7 @@ GmshReader::__proceedData()
             Face({__quadrangles[f][0], __quadrangles[f][1], __quadrangles[f][2], __quadrangles[f][3]},
                  node_number_vector));
           if (i == face_to_id_map.end()) {
-            std::cerr << "face not found!\n";
-            std::terminate();
+            throw NormalError("face not found");
           }
           return i->second;
         }();
@@ -1303,8 +1302,9 @@ GmshReader::__proceedData()
         const unsigned int edge_id = [&] {
           auto i = edge_to_id_map.find(Edge({__edges[e][0], __edges[e][1]}, node_number_vector));
           if (i == edge_to_id_map.end()) {
-            std::cerr << "edge " << __edges[e][0] << " not found!\n";
-            std::terminate();
+            std::stringstream error_msg;
+            error_msg << "edge " << __edges[e][0] << " not found";
+            throw NormalError(error_msg.str());
           }
           return i->second;
         }();
@@ -1461,8 +1461,9 @@ GmshReader::__proceedData()
       const unsigned int edge_id = [&] {
         auto i = face_to_id_map.find(Face({__edges[e][0], __edges[e][1]}, node_number_vector));
         if (i == face_to_id_map.end()) {
-          std::cerr << "face " << __edges[e][0] << " not found!\n";
-          std::terminate();
+          std::stringstream error_msg;
+          error_msg << "face " << __edges[e][0] << " not found";
+          throw NormalError(error_msg.str());
         }
         return i->second;
       }();
diff --git a/src/mesh/ItemValueUtils.hpp b/src/mesh/ItemValueUtils.hpp
index b9570fc9c2bb62d26a8ab2ee1a6b0f9709ba60ac..ca0c38ecdbb4ec48271be555241fbd45f02c66e4 100644
--- a/src/mesh/ItemValueUtils.hpp
+++ b/src/mesh/ItemValueUtils.hpp
@@ -38,8 +38,7 @@ min(const ItemValue<DataType, item_type>& item_value)
       break;
     }
     default: {
-      std::cerr << __FILE__ << ':' << __LINE__ << ": unexpected dimension\n";
-      std::terminate();
+      throw UnexpectedError("unexpected dimension");
     }
     }
   }(*item_value.connectivity_ptr());
@@ -129,8 +128,7 @@ max(const ItemValue<DataType, item_type>& item_value)
       break;
     }
     default: {
-      std::cerr << __FILE__ << ':' << __LINE__ << ": unexpected dimension\n";
-      std::terminate();
+      throw UnexpectedError("unexpected dimension");
     }
     }
   }(*item_value.connectivity_ptr());
@@ -220,8 +218,7 @@ sum(const ItemValue<DataType, item_type>& item_value)
       break;
     }
     default: {
-      std::cerr << __FILE__ << ':' << __LINE__ << ": unexpected dimension\n";
-      std::terminate();
+      throw UnexpectedError("unexpected dimension");
     }
     }
   }(*item_value.connectivity_ptr());
diff --git a/src/mesh/Synchronizer.hpp b/src/mesh/Synchronizer.hpp
index 69b7c8b15bd530941836e366c6db8768b50d93ee..d23634ac3db35a026b50cb4bd0f8c040c80aa4cd 100644
--- a/src/mesh/Synchronizer.hpp
+++ b/src/mesh/Synchronizer.hpp
@@ -3,7 +3,7 @@
 
 #include <mesh/Connectivity.hpp>
 #include <mesh/ItemValue.hpp>
-
+#include <utils/Exceptions.hpp>
 #include <utils/Messenger.hpp>
 
 #include <iostream>
@@ -190,8 +190,7 @@ class Synchronizer
       break;
     }
     default: {
-      std::cerr << __FILE__ << ':' << __LINE__ << ": unexpected dimension\n";
-      std::terminate();
+      throw UnexpectedError("unexpected dimension");
     }
     }
   }