From c64d0a0bdd3abb378d000d9f378f306598e478b9 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 21 Apr 2020 11:37:35 +0200
Subject: [PATCH] Improve exception management.

Remove almost all std::exit calls.

Now, one should use one of the provided exceptions.

- NormalError: an error which is related to the use of the code. It is
not related to a bug

- NotImplementedError: an error related to a functionnality that is
not avaliable yet

- UnexpectedError: an error related to some unexpected state. This
is related to a bug

Closes #13
---
 src/main.cpp                      | 591 +++++++++++++++---------------
 src/mesh/ConnectivityComputer.cpp |   5 +-
 src/mesh/GmshReader.cpp           | 295 ++++++---------
 src/mesh/MeshNodeBoundary.hpp     |  18 +-
 src/output/VTKWriter.hpp          |   7 +-
 src/scheme/AcousticSolver.hpp     |  17 +-
 src/utils/CMakeLists.txt          |   1 +
 src/utils/CastArray.hpp           |   5 +-
 src/utils/Exceptions.cpp          |  28 ++
 src/utils/Exceptions.hpp          |  21 ++
 src/utils/Messenger.cpp           |   3 +-
 src/utils/Messenger.hpp           |   5 +-
 src/utils/Partitioner.cpp         |   5 +-
 src/utils/PugsAssert.hpp          |  29 +-
 src/utils/SignalManager.cpp       |   6 +
 15 files changed, 512 insertions(+), 524 deletions(-)
 create mode 100644 src/utils/Exceptions.cpp
 create mode 100644 src/utils/Exceptions.hpp

diff --git a/src/main.cpp b/src/main.cpp
index 93ff1e39f..346d6ed62 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -10,6 +10,7 @@
 
 #include <VTKWriter.hpp>
 
+#include <Exceptions.hpp>
 #include <Timer.hpp>
 
 #include <TinyMatrix.hpp>
@@ -33,383 +34,387 @@
 int
 main(int argc, char* argv[])
 {
-  std::string filename = initialize(argc, argv);
+  try {
+    std::string filename = initialize(argc, argv);
 
-  std::regex gmsh_regex("(.*).msh");
-  if (not std::regex_match(filename, gmsh_regex)) {
-    parser(filename);
-    return 0;
-  }
+    std::regex gmsh_regex("(.*).msh");
+    if (not std::regex_match(filename, gmsh_regex)) {
+      parser(filename);
+      return 0;
+    }
 
-  std::map<std::string, double> method_cost_map;
+    std::map<std::string, double> method_cost_map;
 
-  SynchronizerManager::create();
+    SynchronizerManager::create();
 
-  if (filename != "") {
-    std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
-    Timer gmsh_timer;
-    gmsh_timer.reset();
-    GmshReader gmsh_reader(filename);
-    method_cost_map["Mesh building"] = gmsh_timer.seconds();
+    if (filename != "") {
+      std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
+      Timer gmsh_timer;
+      gmsh_timer.reset();
+      GmshReader gmsh_reader(filename);
+      method_cost_map["Mesh building"] = gmsh_timer.seconds();
 
-    std::shared_ptr<IMesh> p_mesh = gmsh_reader.mesh();
+      std::shared_ptr<IMesh> p_mesh = gmsh_reader.mesh();
 
-    switch (p_mesh->dimension()) {
-    case 1: {
-      std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
-      std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
-      for (const auto& sym_boundary_name : sym_boundary_name_list) {
-        std::shared_ptr<BoundaryDescriptor> boudary_descriptor =
-          std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
-        SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
-          new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+      switch (p_mesh->dimension()) {
+      case 1: {
+        std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
+        std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
+        for (const auto& sym_boundary_name : sym_boundary_name_list) {
+          std::shared_ptr<BoundaryDescriptor> boudary_descriptor =
+            std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+          SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
+            new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
 
-        bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
-      }
+          bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+        }
 
-      using ConnectivityType = Connectivity1D;
-      using MeshType         = Mesh<ConnectivityType>;
-      using MeshDataType     = MeshData<MeshType>;
-      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
-
-      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
-
-      Timer timer;
-      timer.reset();
-      MeshDataType mesh_data(mesh);
-
-      std::vector<BoundaryConditionHandler> bc_list;
-      {
-        for (const auto& bc_descriptor : bc_descriptor_list) {
-          switch (bc_descriptor->type()) {
-          case BoundaryConditionDescriptor::Type::symmetry: {
-            const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-              dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-            for (size_t i_ref_node_list = 0;
-                 i_ref_node_list < mesh.connectivity().numberOfRefItemList<ItemType::node>(); ++i_ref_node_list) {
-              const RefNodeList& ref_node_list = mesh.connectivity().refItemList<ItemType::node>(i_ref_node_list);
-              const RefId& ref                 = ref_node_list.refId();
-              if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-                SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
-                  new SymmetryBoundaryCondition<MeshType::Dimension>(
-                    MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_node_list));
-                std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
-                bc_list.push_back(BoundaryConditionHandler(bc));
+        using ConnectivityType = Connectivity1D;
+        using MeshType         = Mesh<ConnectivityType>;
+        using MeshDataType     = MeshData<MeshType>;
+        using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+        const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+        Timer timer;
+        timer.reset();
+        MeshDataType mesh_data(mesh);
+
+        std::vector<BoundaryConditionHandler> bc_list;
+        {
+          for (const auto& bc_descriptor : bc_descriptor_list) {
+            switch (bc_descriptor->type()) {
+            case BoundaryConditionDescriptor::Type::symmetry: {
+              const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+                dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+              for (size_t i_ref_node_list = 0;
+                   i_ref_node_list < mesh.connectivity().numberOfRefItemList<ItemType::node>(); ++i_ref_node_list) {
+                const RefNodeList& ref_node_list = mesh.connectivity().refItemList<ItemType::node>(i_ref_node_list);
+                const RefId& ref                 = ref_node_list.refId();
+                if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+                  SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
+                    new SymmetryBoundaryCondition<MeshType::Dimension>(
+                      MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_node_list));
+                  std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
+                  bc_list.push_back(BoundaryConditionHandler(bc));
+                }
               }
+              break;
+            }
+            default: {
+              throw UnexpectedError("Unknown BCDescription\n");
+            }
             }
-            break;
-          }
-          default: {
-            std::cerr << "Unknown BCDescription\n";
-            std::exit(1);
-          }
           }
         }
-      }
 
-      UnknownsType unknowns(mesh_data);
+        UnknownsType unknowns(mesh_data);
+
+        unknowns.initializeSod();
 
-      unknowns.initializeSod();
+        AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
 
-      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+        using Rd = TinyVector<MeshType::Dimension>;
 
-      using Rd = TinyVector<MeshType::Dimension>;
+        const CellValue<const double>& Vj = mesh_data.Vj();
 
-      const CellValue<const double>& Vj = mesh_data.Vj();
+        const double tmax = 0.2;
+        double t          = 0;
 
-      const double tmax = 0.2;
-      double t          = 0;
+        int itermax   = std::numeric_limits<int>::max();
+        int iteration = 0;
 
-      int itermax   = std::numeric_limits<int>::max();
-      int iteration = 0;
+        CellValue<double>& rhoj   = unknowns.rhoj();
+        CellValue<double>& ej     = unknowns.ej();
+        CellValue<double>& pj     = unknowns.pj();
+        CellValue<double>& gammaj = unknowns.gammaj();
+        CellValue<double>& cj     = unknowns.cj();
 
-      CellValue<double>& rhoj   = unknowns.rhoj();
-      CellValue<double>& ej     = unknowns.ej();
-      CellValue<double>& pj     = unknowns.pj();
-      CellValue<double>& gammaj = unknowns.gammaj();
-      CellValue<double>& cj     = unknowns.cj();
+        BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
-      BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+        VTKWriter vtk_writer("mesh", 0.01);
 
-      VTKWriter vtk_writer("mesh", 0.01);
+        while ((t < tmax) and (iteration < itermax)) {
+          vtk_writer.write(mesh,
+                           {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                            NamedItemValue{"coords", mesh.xr()},
+                            NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+                            NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                           t);
+          double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
+          if (t + dt > tmax) {
+            dt = tmax - t;
+          }
+          acoustic_solver.computeNextStep(t, dt, unknowns);
+
+          block_eos.updatePandCFromRhoE();
 
-      while ((t < tmax) and (iteration < itermax)) {
+          t += dt;
+          ++iteration;
+        }
         vtk_writer.write(mesh,
                          {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                           NamedItemValue{"coords", mesh.xr()},
                           NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                           NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                         t);
-        double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
-        if (t + dt > tmax) {
-          dt = tmax - t;
-        }
-        acoustic_solver.computeNextStep(t, dt, unknowns);
+                         t, true);   // forces last output
 
-        block_eos.updatePandCFromRhoE();
+        std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
+                  << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
-        t += dt;
-        ++iteration;
-      }
-      vtk_writer.write(mesh,
-                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", mesh.xr()},
-                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                       t, true);   // forces last output
-
-      std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
-                << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
-      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
-
-      {   // gnuplot output for density
-        const CellValue<const Rd>& xj       = mesh_data.xj();
-        const CellValue<const double>& rhoj = unknowns.rhoj();
-        std::ofstream fout("rho");
-        for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
-          fout << xj[j][0] << ' ' << rhoj[j] << '\n';
+        method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+
+        {   // gnuplot output for density
+          const CellValue<const Rd>& xj       = mesh_data.xj();
+          const CellValue<const double>& rhoj = unknowns.rhoj();
+          std::ofstream fout("rho");
+          for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+            fout << xj[j][0] << ' ' << rhoj[j] << '\n';
+          }
         }
-      }
 
-      break;
-    }
-    case 2: {
-      // test case boundary condition description
-      std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX"};
-      std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
-      for (const auto& sym_boundary_name : sym_boundary_name_list) {
-        std::shared_ptr<BoundaryDescriptor> boudary_descriptor =
-          std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
-        SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
-          new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
-        bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+        break;
       }
+      case 2: {
+        // test case boundary condition description
+        std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX"};
+        std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
+        for (const auto& sym_boundary_name : sym_boundary_name_list) {
+          std::shared_ptr<BoundaryDescriptor> boudary_descriptor =
+            std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+          SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
+            new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+          bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+        }
 
-      using ConnectivityType = Connectivity2D;
-      using MeshType         = Mesh<ConnectivityType>;
-      using MeshDataType     = MeshData<MeshType>;
-      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
-
-      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
-
-      Timer timer;
-      timer.reset();
-      MeshDataType mesh_data(mesh);
-
-      std::vector<BoundaryConditionHandler> bc_list;
-      {
-        for (const auto& bc_descriptor : bc_descriptor_list) {
-          switch (bc_descriptor->type()) {
-          case BoundaryConditionDescriptor::Type::symmetry: {
-            const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-              dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-            for (size_t i_ref_face_list = 0;
-                 i_ref_face_list < mesh.connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
-              const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
-              const RefId& ref                 = ref_face_list.refId();
-              if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-                SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
-                  new SymmetryBoundaryCondition<MeshType::Dimension>(
-                    MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
-                std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
-                bc_list.push_back(BoundaryConditionHandler(bc));
+        using ConnectivityType = Connectivity2D;
+        using MeshType         = Mesh<ConnectivityType>;
+        using MeshDataType     = MeshData<MeshType>;
+        using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+        const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+        Timer timer;
+        timer.reset();
+        MeshDataType mesh_data(mesh);
+
+        std::vector<BoundaryConditionHandler> bc_list;
+        {
+          for (const auto& bc_descriptor : bc_descriptor_list) {
+            switch (bc_descriptor->type()) {
+            case BoundaryConditionDescriptor::Type::symmetry: {
+              const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+                dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+              for (size_t i_ref_face_list = 0;
+                   i_ref_face_list < mesh.connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
+                const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
+                const RefId& ref                 = ref_face_list.refId();
+                if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+                  SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
+                    new SymmetryBoundaryCondition<MeshType::Dimension>(
+                      MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
+                  std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
+                  bc_list.push_back(BoundaryConditionHandler(bc));
+                }
               }
+              break;
+            }
+            default: {
+              throw UnexpectedError("Unknown BCDescription\n");
+            }
             }
-            break;
-          }
-          default: {
-            std::cerr << "Unknown BCDescription\n";
-            std::exit(1);
-          }
           }
         }
-      }
 
-      UnknownsType unknowns(mesh_data);
+        UnknownsType unknowns(mesh_data);
+
+        unknowns.initializeSod();
 
-      unknowns.initializeSod();
+        AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
 
-      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+        const CellValue<const double>& Vj = mesh_data.Vj();
 
-      const CellValue<const double>& Vj = mesh_data.Vj();
+        const double tmax = 0.2;
+        double t          = 0;
 
-      const double tmax = 0.2;
-      double t          = 0;
+        int itermax   = std::numeric_limits<int>::max();
+        int iteration = 0;
 
-      int itermax   = std::numeric_limits<int>::max();
-      int iteration = 0;
+        CellValue<double>& rhoj   = unknowns.rhoj();
+        CellValue<double>& ej     = unknowns.ej();
+        CellValue<double>& pj     = unknowns.pj();
+        CellValue<double>& gammaj = unknowns.gammaj();
+        CellValue<double>& cj     = unknowns.cj();
 
-      CellValue<double>& rhoj   = unknowns.rhoj();
-      CellValue<double>& ej     = unknowns.ej();
-      CellValue<double>& pj     = unknowns.pj();
-      CellValue<double>& gammaj = unknowns.gammaj();
-      CellValue<double>& cj     = unknowns.cj();
+        BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
-      BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+        VTKWriter vtk_writer("mesh", 0.01);
 
-      VTKWriter vtk_writer("mesh", 0.01);
+        while ((t < tmax) and (iteration < itermax)) {
+          vtk_writer.write(mesh,
+                           {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                            NamedItemValue{"coords", mesh.xr()},
+                            NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+                            NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                           t);
+          double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
+          if (t + dt > tmax) {
+            dt = tmax - t;
+          }
+          acoustic_solver.computeNextStep(t, dt, unknowns);
+
+          block_eos.updatePandCFromRhoE();
 
-      while ((t < tmax) and (iteration < itermax)) {
+          t += dt;
+          ++iteration;
+        }
         vtk_writer.write(mesh,
                          {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                           NamedItemValue{"coords", mesh.xr()},
                           NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                           NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                         t);
-        double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
-        if (t + dt > tmax) {
-          dt = tmax - t;
-        }
-        acoustic_solver.computeNextStep(t, dt, unknowns);
+                         t, true);   // forces last output
 
-        block_eos.updatePandCFromRhoE();
+        std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
+                  << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
-        t += dt;
-        ++iteration;
-      }
-      vtk_writer.write(mesh,
-                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", mesh.xr()},
-                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                       t, true);   // forces last output
-
-      std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
-                << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
-      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
-      break;
-    }
-    case 3: {
-      std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
-      std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
-      for (const auto& sym_boundary_name : sym_boundary_name_list) {
-        std::shared_ptr<BoundaryDescriptor> boudary_descriptor =
-          std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
-        SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
-          new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
-        bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+        method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+        break;
       }
+      case 3: {
+        std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
+        std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
+        for (const auto& sym_boundary_name : sym_boundary_name_list) {
+          std::shared_ptr<BoundaryDescriptor> boudary_descriptor =
+            std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+          SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
+            new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+          bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+        }
 
-      using ConnectivityType = Connectivity3D;
-      using MeshType         = Mesh<ConnectivityType>;
-      using MeshDataType     = MeshData<MeshType>;
-      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
-
-      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
-
-      Timer timer;
-      timer.reset();
-      MeshDataType mesh_data(mesh);
-
-      std::vector<BoundaryConditionHandler> bc_list;
-      {
-        for (const auto& bc_descriptor : bc_descriptor_list) {
-          switch (bc_descriptor->type()) {
-          case BoundaryConditionDescriptor::Type::symmetry: {
-            const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-              dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-            for (size_t i_ref_face_list = 0;
-                 i_ref_face_list < mesh.connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
-              const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
-              const RefId& ref                 = ref_face_list.refId();
-              if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-                SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
-                  new SymmetryBoundaryCondition<MeshType::Dimension>(
-                    MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
-                std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
-                bc_list.push_back(BoundaryConditionHandler(bc));
+        using ConnectivityType = Connectivity3D;
+        using MeshType         = Mesh<ConnectivityType>;
+        using MeshDataType     = MeshData<MeshType>;
+        using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+        const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+        Timer timer;
+        timer.reset();
+        MeshDataType mesh_data(mesh);
+
+        std::vector<BoundaryConditionHandler> bc_list;
+        {
+          for (const auto& bc_descriptor : bc_descriptor_list) {
+            switch (bc_descriptor->type()) {
+            case BoundaryConditionDescriptor::Type::symmetry: {
+              const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+                dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+              for (size_t i_ref_face_list = 0;
+                   i_ref_face_list < mesh.connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
+                const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
+                const RefId& ref                 = ref_face_list.refId();
+                if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+                  SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
+                    new SymmetryBoundaryCondition<MeshType::Dimension>(
+                      MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
+                  std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
+                  bc_list.push_back(BoundaryConditionHandler(bc));
+                }
               }
+              break;
+            }
+            default: {
+              throw UnexpectedError("Unknown BCDescription\n");
+            }
             }
-            break;
-          }
-          default: {
-            std::cerr << "Unknown BCDescription\n";
-            std::exit(1);
-          }
           }
         }
-      }
 
-      UnknownsType unknowns(mesh_data);
+        UnknownsType unknowns(mesh_data);
+
+        unknowns.initializeSod();
 
-      unknowns.initializeSod();
+        AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
 
-      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+        const CellValue<const double>& Vj = mesh_data.Vj();
 
-      const CellValue<const double>& Vj = mesh_data.Vj();
+        const double tmax = 0.2;
+        double t          = 0;
 
-      const double tmax = 0.2;
-      double t          = 0;
+        int itermax   = std::numeric_limits<int>::max();
+        int iteration = 0;
 
-      int itermax   = std::numeric_limits<int>::max();
-      int iteration = 0;
+        CellValue<double>& rhoj   = unknowns.rhoj();
+        CellValue<double>& ej     = unknowns.ej();
+        CellValue<double>& pj     = unknowns.pj();
+        CellValue<double>& gammaj = unknowns.gammaj();
+        CellValue<double>& cj     = unknowns.cj();
 
-      CellValue<double>& rhoj   = unknowns.rhoj();
-      CellValue<double>& ej     = unknowns.ej();
-      CellValue<double>& pj     = unknowns.pj();
-      CellValue<double>& gammaj = unknowns.gammaj();
-      CellValue<double>& cj     = unknowns.cj();
+        BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
-      BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+        VTKWriter vtk_writer("mesh", 0.01);
 
-      VTKWriter vtk_writer("mesh", 0.01);
+        while ((t < tmax) and (iteration < itermax)) {
+          vtk_writer.write(mesh,
+                           {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                            NamedItemValue{"coords", mesh.xr()},
+                            NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+                            NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                           t);
+          double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
+          if (t + dt > tmax) {
+            dt = tmax - t;
+          }
+          acoustic_solver.computeNextStep(t, dt, unknowns);
+          block_eos.updatePandCFromRhoE();
 
-      while ((t < tmax) and (iteration < itermax)) {
+          t += dt;
+          ++iteration;
+        }
         vtk_writer.write(mesh,
                          {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                           NamedItemValue{"coords", mesh.xr()},
                           NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                           NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                         t);
-        double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
-        if (t + dt > tmax) {
-          dt = tmax - t;
-        }
-        acoustic_solver.computeNextStep(t, dt, unknowns);
-        block_eos.updatePandCFromRhoE();
+                         t, true);   // forces last output
+
+        std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
+                  << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
-        t += dt;
-        ++iteration;
+        method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+        break;
+      }
       }
-      vtk_writer.write(mesh,
-                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", mesh.xr()},
-                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                       t, true);   // forces last output
-
-      std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
-                << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
-      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
-      break;
-    }
-    }
 
-    std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':'
-              << __LINE__ << ")\n";
+      std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':'
+                << __LINE__ << ")\n";
 
-  } else {
-    std::cerr << "Connectivity1D defined by number of nodes no more implemented\n";
-    std::exit(0);
-  }
+    } else {
+      throw NormalError("Connectivity1D defined by number of nodes no more implemented\n");
+    }
 
-  SynchronizerManager::destroy();
+    SynchronizerManager::destroy();
 
-  finalize();
+    finalize();
 
-  std::string::size_type size = 0;
-  for (const auto& method_cost : method_cost_map) {
-    size = std::max(size, method_cost.first.size());
-  }
+    std::string::size_type size = 0;
+    for (const auto& method_cost : method_cost_map) {
+      size = std::max(size, method_cost.first.size());
+    }
 
-  for (const auto& method_cost : method_cost_map) {
-    std::cout << "* [" << rang::fgB::cyan << std::setw(size) << std::left << method_cost.first << rang::fg::reset
-              << "] Execution time: " << rang::style::bold << method_cost.second << rang::style::reset << '\n';
+    for (const auto& method_cost : method_cost_map) {
+      std::cout << "* [" << rang::fgB::cyan << std::setw(size) << std::left << method_cost.first << rang::fg::reset
+                << "] Execution time: " << rang::style::bold << method_cost.second << rang::style::reset << '\n';
+    }
+  }
+  catch (const NormalError& e) {
+    // Each failing process must write
+    std::cerr.setstate(std::ios::goodbit);
+    std::cerr << e.what() << '\n';
+    return 1;
   }
 
   return 0;
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index b4f24fb69..8e1d5d526 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -2,6 +2,8 @@
 
 #include <ConnectivityComputer.hpp>
 
+#include <Exceptions.hpp>
+
 #include <iostream>
 #include <map>
 
@@ -56,8 +58,7 @@ ConnectivityComputer::_computeInverse(const ConnectivityMatrix& item_to_child_ma
     size_t i = 0;
     for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
       if (child_item_id != i) {
-        std::cerr << "sparse item numerotation NIY\n";
-        std::exit(0);
+        throw NotImplementedError("sparse item numerotation NIY\n");
       }
       ++i;
     }
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 7ba1dbfc1..25a224caf 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -20,8 +20,11 @@
 
 #include <ConnectivityDispatcher.hpp>
 
+#include <Exceptions.hpp>
+
 #include <iomanip>
 #include <iostream>
+#include <sstream>
 
 #include <map>
 #include <regex>
@@ -342,160 +345,115 @@ class ConnectivityFace<3>
 GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
 {
   if (parallel::rank() == 0) {
-    try {
-      m_fin.open(m_filename);
-      if (not m_fin) {
-        std::cerr << "cannot read file '" << m_filename << "'\n";
-        std::exit(0);
+    m_fin.open(m_filename);
+    if (not m_fin) {
+      std::ostringstream os;
+      os << "cannot read file '" << m_filename << "'";
+      throw NormalError(os.str());
+    }
+
+    // gmsh 2.2 format keywords
+    __keywordList["$MeshFormat"]       = MESHFORMAT;
+    __keywordList["$EndMeshFormat"]    = ENDMESHFORMAT;
+    __keywordList["$Nodes"]            = NODES;
+    __keywordList["$EndNodes"]         = ENDNODES;
+    __keywordList["$Elements"]         = ELEMENTS;
+    __keywordList["$EndElements"]      = ENDELEMENTS;
+    __keywordList["$PhysicalNames"]    = PHYSICALNAMES;
+    __keywordList["$EndPhysicalNames"] = ENDPHYSICALNAMES;
+
+    __numberOfPrimitiveNodes.resize(16);
+    __numberOfPrimitiveNodes[0]  = 2;    // edge
+    __numberOfPrimitiveNodes[1]  = 3;    // triangle
+    __numberOfPrimitiveNodes[2]  = 4;    // quadrangle
+    __numberOfPrimitiveNodes[3]  = 4;    // Tetrahedron
+    __numberOfPrimitiveNodes[4]  = 8;    // Hexaredron
+    __numberOfPrimitiveNodes[5]  = 6;    // Prism
+    __numberOfPrimitiveNodes[6]  = 5;    // Pyramid
+    __numberOfPrimitiveNodes[7]  = 3;    // second order edge
+    __numberOfPrimitiveNodes[8]  = 6;    // second order triangle
+    __numberOfPrimitiveNodes[9]  = 9;    // second order quadrangle
+    __numberOfPrimitiveNodes[10] = 10;   // second order tetrahedron
+    __numberOfPrimitiveNodes[11] = 27;   // second order hexahedron
+    __numberOfPrimitiveNodes[12] = 18;   // second order prism
+    __numberOfPrimitiveNodes[13] = 14;   // second order pyramid
+    __numberOfPrimitiveNodes[14] = 1;    // point
+
+    __primitivesNames[0]      = "edges";
+    __supportedPrimitives[0]  = true;
+    __primitivesNames[1]      = "triangles";
+    __supportedPrimitives[1]  = true;
+    __primitivesNames[2]      = "quadrangles";
+    __supportedPrimitives[2]  = true;
+    __primitivesNames[3]      = "tetrahedra";
+    __supportedPrimitives[3]  = true;
+    __primitivesNames[4]      = "hexahedra";
+    __supportedPrimitives[4]  = true;
+    __primitivesNames[5]      = "prisms";
+    __supportedPrimitives[5]  = false;
+    __primitivesNames[6]      = "pyramids";
+    __supportedPrimitives[6]  = false;
+    __primitivesNames[7]      = "second order edges";
+    __supportedPrimitives[7]  = false;
+    __primitivesNames[8]      = "second order triangles";
+    __supportedPrimitives[8]  = false;
+    __primitivesNames[9]      = "second order quadrangles";
+    __supportedPrimitives[9]  = false;
+    __primitivesNames[10]     = "second order tetrahedra";
+    __supportedPrimitives[10] = false;
+    __primitivesNames[11]     = "second order hexahedra";
+    __supportedPrimitives[11] = false;
+    __primitivesNames[12]     = "second order prisms";
+    __supportedPrimitives[12] = false;
+    __primitivesNames[13]     = "second order pyramids";
+    __supportedPrimitives[13] = false;
+    __primitivesNames[14]     = "point";
+    __supportedPrimitives[14] = true;
+
+    std::cout << "Reading file '" << m_filename << "'\n";
+
+    // Getting vertices list
+    GmshReader::Keyword kw = this->__nextKeyword();
+    switch (kw.second) {
+    // case NOD: {
+    //   this->__readGmsh1();
+    //   break;
+    // }
+    case MESHFORMAT: {
+      double fileVersion = this->_getReal();
+      if (fileVersion != 2.2) {
+        throw NormalError("Cannot read Gmsh format '" + std::to_string(fileVersion) + "'");
+      }
+      int fileType = this->_getInteger();
+      __binary     = (fileType == 1);
+      if ((fileType < 0) or (fileType > 1)) {
+        throw NormalError("Cannot read Gmsh file type '" + std::to_string(fileType) + "'");
       }
 
-      // gmsh 2.2 format keywords
-      __keywordList["$MeshFormat"]       = MESHFORMAT;
-      __keywordList["$EndMeshFormat"]    = ENDMESHFORMAT;
-      __keywordList["$Nodes"]            = NODES;
-      __keywordList["$EndNodes"]         = ENDNODES;
-      __keywordList["$Elements"]         = ELEMENTS;
-      __keywordList["$EndElements"]      = ENDELEMENTS;
-      __keywordList["$PhysicalNames"]    = PHYSICALNAMES;
-      __keywordList["$EndPhysicalNames"] = ENDPHYSICALNAMES;
-
-      __numberOfPrimitiveNodes.resize(16);
-      __numberOfPrimitiveNodes[0]  = 2;    // edge
-      __numberOfPrimitiveNodes[1]  = 3;    // triangle
-      __numberOfPrimitiveNodes[2]  = 4;    // quadrangle
-      __numberOfPrimitiveNodes[3]  = 4;    // Tetrahedron
-      __numberOfPrimitiveNodes[4]  = 8;    // Hexaredron
-      __numberOfPrimitiveNodes[5]  = 6;    // Prism
-      __numberOfPrimitiveNodes[6]  = 5;    // Pyramid
-      __numberOfPrimitiveNodes[7]  = 3;    // second order edge
-      __numberOfPrimitiveNodes[8]  = 6;    // second order triangle
-      __numberOfPrimitiveNodes[9]  = 9;    // second order quadrangle
-      __numberOfPrimitiveNodes[10] = 10;   // second order tetrahedron
-      __numberOfPrimitiveNodes[11] = 27;   // second order hexahedron
-      __numberOfPrimitiveNodes[12] = 18;   // second order prism
-      __numberOfPrimitiveNodes[13] = 14;   // second order pyramid
-      __numberOfPrimitiveNodes[14] = 1;    // point
-
-      __primitivesNames[0]      = "edges";
-      __supportedPrimitives[0]  = true;
-      __primitivesNames[1]      = "triangles";
-      __supportedPrimitives[1]  = true;
-      __primitivesNames[2]      = "quadrangles";
-      __supportedPrimitives[2]  = true;
-      __primitivesNames[3]      = "tetrahedra";
-      __supportedPrimitives[3]  = true;
-      __primitivesNames[4]      = "hexahedra";
-      __supportedPrimitives[4]  = true;
-      __primitivesNames[5]      = "prisms";
-      __supportedPrimitives[5]  = false;
-      __primitivesNames[6]      = "pyramids";
-      __supportedPrimitives[6]  = false;
-      __primitivesNames[7]      = "second order edges";
-      __supportedPrimitives[7]  = false;
-      __primitivesNames[8]      = "second order triangles";
-      __supportedPrimitives[8]  = false;
-      __primitivesNames[9]      = "second order quadrangles";
-      __supportedPrimitives[9]  = false;
-      __primitivesNames[10]     = "second order tetrahedra";
-      __supportedPrimitives[10] = false;
-      __primitivesNames[11]     = "second order hexahedra";
-      __supportedPrimitives[11] = false;
-      __primitivesNames[12]     = "second order prisms";
-      __supportedPrimitives[12] = false;
-      __primitivesNames[13]     = "second order pyramids";
-      __supportedPrimitives[13] = false;
-      __primitivesNames[14]     = "point";
-      __supportedPrimitives[14] = true;
-
-      std::cout << "Reading file '" << m_filename << "'\n";
-
-      // Getting vertices list
-      GmshReader::Keyword kw = this->__nextKeyword();
-      switch (kw.second) {
-      // case NOD: {
-      //   this->__readGmsh1();
-      //   break;
-      // }
-      case MESHFORMAT: {
-        double fileVersion = this->_getReal();
-        if (fileVersion != 2.2) {
-          throw ErrorHandler(__FILE__, __LINE__, "Cannot read Gmsh format '" + std::to_string(fileVersion) + "'",
-                             ErrorHandler::normal);
-        }
-        int fileType = this->_getInteger();
-        __binary     = (fileType == 1);
-        if ((fileType < 0) or (fileType > 1)) {
-          throw ErrorHandler(__FILE__, __LINE__, "Cannot read Gmsh file type '" + std::to_string(fileType) + "'",
-                             ErrorHandler::normal);
-        }
-
-        int dataSize = this->_getInteger();
-        if (dataSize != sizeof(double)) {
-          throw ErrorHandler(__FILE__, __LINE__, "Data size not supported '" + std::to_string(dataSize) + "'",
-                             ErrorHandler::normal);
-        }
-
-        if (__binary) {
-          //       int one=0;
-
-          //       fseek(__ifh,1L,SEEK_CUR);
-          //       fread(reinterpret_cast<char*>(&one),sizeof(int),1,__ifh);
-
-          //       if (one == 1) {
-          // 	__convertEndian = false;
-          //       } else {
-          // 	invertEndianess(one);
-
-          // 	if (one == 1) {
-          // 	  __convertEndian = true;
-          // 	} else {
-          // 	  throw ErrorHandler(__FILE__,__LINE__,
-          // 	  		     "Cannot determine endianess",
-          // 	  		     ErrorHandler::normal);
-          // 	}
-          // }
-
-          //       std::cout << "- Binary format: ";
-          // #ifdef    WORDS_BIGENDIAN
-          //       if (not __convertEndian) {
-          // 	std::cout << "big endian\n";
-          //       } else {
-          // 	std::cout << "little endian\n";
-          //       }
-          // #else  // WORDS_BIGENDIAN
-          //       if (not __convertEndian) {
-          // 	std::cout << "little endian\n";
-          //       } else {
-          // 	std::cout << "big endian\n";
-          //       }
-          // #endif // WORDS_BIGENDIAN
-        }
-
-        kw = this->__nextKeyword();
-        if (kw.second != ENDMESHFORMAT) {
-          throw ErrorHandler(__FILE__, __LINE__,
-                             "reading file '" + m_filename + "': expecting $EndMeshFormat, '" + kw.first +
-                               "' was found",
-                             ErrorHandler::normal);
-        }
-
-        this->__readGmshFormat2_2();
-
-        break;
+      int dataSize = this->_getInteger();
+      if (dataSize != sizeof(double)) {
+        throw NormalError("Data size not supported '" + std::to_string(dataSize) + "'");
       }
-      default: {
-        throw ErrorHandler(__FILE__, __LINE__, "cannot determine format version of '" + m_filename + "'",
-                           ErrorHandler::normal);
+
+      if (__binary) {
+        throw NotImplementedError("cannot read binary files");
       }
+
+      kw = this->__nextKeyword();
+      if (kw.second != ENDMESHFORMAT) {
+        throw NormalError("reading file '" + m_filename + "': expecting $EndMeshFormat, '" + kw.first + "' was found");
       }
 
-      this->__proceedData();
-      // this->__createMesh();
+      this->__readGmshFormat2_2();
+
+      break;
     }
-    catch (const ErrorHandler& e) {
-      e.writeErrorMessage();
-      std::exit(0);
+    default: {
+      throw NormalError("cannot determine format version of '" + m_filename + "'");
     }
+    }
+
+    this->__proceedData();
   }
   std::cout << std::flush;
   if (parallel::size() > 1) {
@@ -516,8 +474,7 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
         }
       }
       if (dimension_set.size() != 1) {
-        std::cerr << "error dimensions of read mesh parts differ!\n";
-        std::exit(1);
+        throw NormalError("dimensions of read mesh parts differ!\n");
       }
 
       return *begin(dimension_set);
@@ -536,8 +493,7 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
       break;
     }
     default: {
-      std::cerr << "unexpected dimension " << mesh_dimension << '\n';
-      std::exit(1);
+      throw NormalError("invalid mesh dimension" + std::to_string((mesh_dimension)));
     }
     }
   }
@@ -549,8 +505,7 @@ GmshReader::__readVertices()
   const int numberOfVerices = this->_getInteger();
   std::cout << "- Number Of Vertices: " << numberOfVerices << '\n';
   if (numberOfVerices < 0) {
-    throw ErrorHandler(__FILE__, __LINE__, "reading file '" + this->m_filename + "': number of vertices is negative",
-                       ErrorHandler::normal);
+    throw NormalError("reading file '" + this->m_filename + "': number of vertices is negative");
   }
 
   __verticesNumbers.resize(numberOfVerices);
@@ -565,26 +520,7 @@ GmshReader::__readVertices()
       __vertices[i]        = TinyVector<3, double>(x, y, z);
     }
   } else {
-    // fseek(m_fin.file()__ifh,1L,SEEK_CUR);
-    // for (int i=0; i<numberOfVerices; ++i) {
-    //   int number = 0;
-    //   fread(reinterpret_cast<char*>(&number),sizeof(int),1,__ifh);
-    //   __verticesNumbers[i] = number;
-    //   TinyVector<3,double> x;
-    //   fread(reinterpret_cast<char*>(&(x[0])),sizeof(double),3,__ifh);
-    //   for (size_t j=0; j<3; ++j) {
-    // 	__vertices[i][j] = x[j];
-    //   }
-    // }
-
-    // if (__convertEndian) {
-    //   for (int i=0; i<numberOfVerices; ++i) {
-    // 	invertEndianess(__verticesNumbers[i]);
-    // 	for (size_t j=0; j<3; ++j) {
-    // 	  invertEndianess(vertices[i][j]);
-    // 	}
-    //   }
-    // }
+    throw NotImplementedError("cannot read binary format");
   }
 }
 
@@ -754,11 +690,13 @@ GmshReader::__readPhysicalNames2_2()
     physical_name = std::regex_replace(physical_name, std::regex("(\")"), "");
 
     PhysicalRefId physical_ref_id(physical_dimension, RefId(physical_number, physical_name));
-    auto searched_physical_ref_id = m_physical_ref_map.find(physical_number);
-    if (searched_physical_ref_id != m_physical_ref_map.end()) {
-      std::cerr << "Physical reference id '" << physical_ref_id << "' already defined as '"
-                << searched_physical_ref_id->second << "'!";
-      std::exit(0);
+
+    if (auto i_searched_physical_ref_id = m_physical_ref_map.find(physical_number);
+        i_searched_physical_ref_id != m_physical_ref_map.end()) {
+      std::stringstream os;
+      os << "Physical reference id '" << physical_ref_id << "' already defined as '"
+         << i_searched_physical_ref_id->second << "'!";
+      throw NormalError(os.str());
     }
     m_physical_ref_map[physical_number] = physical_ref_id;
   }
@@ -1778,8 +1716,7 @@ GmshReader::__proceedData()
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
 
   } else {
-    std::cerr << "*** using a dimension 0 mesh is forbidden!\n";
-    std::exit(0);
+    throw NormalError("using a dimension 0 mesh is forbidden");
   }
 }
 
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 5c5de47d4..fdb771eb5 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -13,6 +13,7 @@
 #include <IConnectivity.hpp>
 #include <iostream>
 
+#include <Exceptions.hpp>
 #include <Messenger.hpp>
 
 #include <iostream>
@@ -44,8 +45,7 @@ class MeshNodeBoundary
       face_list.size(), PUGS_LAMBDA(const int& l) {
         const auto& face_cells = face_to_cell_matrix[face_list[l]];
         if (face_cells.size() > 1) {
-          std::cerr << "internal faces cannot be used to define mesh boundaries\n";
-          std::exit(1);
+          throw NormalError("internal faces cannot be used to define mesh boundaries");
         }
       });
 
@@ -155,8 +155,7 @@ MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(const TinyVector<2, double>& norma
     m_node_list.size(), PUGS_LAMBDA(const size_t& r) {
       const R2& x = xr[m_node_list[r]];
       if ((x - origin, normal) > 1E-13 * length) {
-        std::cerr << "this FlatBoundary is not flat!\n";
-        std::exit(1);
+        throw NormalError("this FlatBoundary is not flat!");
       }
     });
 }
@@ -179,8 +178,7 @@ MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
   }();
 
   if (number_of_bc_nodes != 1) {
-    std::cerr << "Node boundaries in 1D require to have exactly 1 node\n";
-    std::exit(1);
+    throw NormalError("Node boundaries in 1D require to have exactly 1 node");
   }
 
   return R{1};
@@ -226,8 +224,9 @@ MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
   }
 
   if (xmin == xmax) {
-    std::cerr << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";
-    std::exit(1);
+    std::stringstream os;
+    os << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";
+    throw NormalError(os.str());
   }
 
   R2 dx = xmax - xmin;
@@ -351,8 +350,7 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
   }
 
   if (normal_l2 == 0) {
-    std::cerr << "not able to compute normal!\n";
-    std::exit(1);
+    throw NormalError("cannot to compute normal!");
   }
 
   normal *= 1. / sqrt(normal_l2);
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index acd10acd5..c48da737a 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -8,6 +8,8 @@
 #include <Messenger.hpp>
 #include <OutputNamedItemValueSet.hpp>
 
+#include <Exceptions.hpp>
+
 #include <fstream>
 #include <iomanip>
 #include <iostream>
@@ -347,8 +349,9 @@ class VTKWriter
               break;
             }
             default: {
-              std::cerr << __FILE__ << ':' << __LINE__ << ": unknown cell type\n";
-              std::exit(1);
+              std::ostringstream os;
+              os << __FILE__ << ':' << __LINE__ << ": unknown cell type";
+              throw UnexpectedError(os.str());
             }
             }
           });
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 24c2e3651..3fb130013 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -19,6 +19,8 @@
 #include <Messenger.hpp>
 #include <SubItemValuePerItem.hpp>
 
+#include <Exceptions.hpp>
+
 #include <iostream>
 
 template <typename MeshData>
@@ -120,22 +122,13 @@ class AcousticSolver
     for (const auto& handler : m_boundary_condition_list) {
       switch (handler.boundaryCondition().type()) {
       case BoundaryCondition::normal_velocity: {
-        std::cerr << __FILE__ << ':' << __LINE__ << ": normal_velocity BC NIY\n";
-        std::exit(0);
-        break;
+        throw NotImplementedError("normal_velocity BC");
       }
       case BoundaryCondition::velocity: {
-        std::cerr << __FILE__ << ':' << __LINE__ << ": velocity BC NIY\n";
-        std::exit(0);
-        break;
+        throw NotImplementedError("velocity BC");
       }
       case BoundaryCondition::pressure: {
-        // const PressureBoundaryCondition& pressure_bc
-        //   = dynamic_cast<const
-        //   PressureBoundaryCondition&>(handler.boundaryCondition());
-        std::cerr << __FILE__ << ':' << __LINE__ << ": pressure BC NIY\n";
-        std::exit(0);
-        break;
+        throw NotImplementedError("pressure BC");
       }
       case BoundaryCondition::symmetry: {
         const SymmetryBoundaryCondition<Dimension>& symmetry_bc =
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index df37bba9e..9c4e30584 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -8,6 +8,7 @@ add_library(
   BuildInfo.cpp
   BacktraceManager.cpp
   ConsoleManager.cpp
+  Exceptions.cpp
   FPEManager.cpp
   Messenger.cpp
   Partitioner.cpp
diff --git a/src/utils/CastArray.hpp b/src/utils/CastArray.hpp
index 5202d8862..b7587abcb 100644
--- a/src/utils/CastArray.hpp
+++ b/src/utils/CastArray.hpp
@@ -6,6 +6,8 @@
 
 #include <iostream>
 
+#include <Exceptions.hpp>
+
 template <typename DataType, typename CastDataType>
 class CastArray
 {
@@ -54,8 +56,7 @@ class CastArray
                   "CastArray cannot remove const attribute");
 
     if (sizeof(DataType) * array.size() % sizeof(CastDataType)) {
-      std::cerr << "cannot cast array to the chosen data type\n";
-      std::exit(1);
+      throw UnexpectedError("cannot cast array to the chosen data type");
     }
   }
 
diff --git a/src/utils/Exceptions.cpp b/src/utils/Exceptions.cpp
new file mode 100644
index 000000000..fa3220716
--- /dev/null
+++ b/src/utils/Exceptions.cpp
@@ -0,0 +1,28 @@
+#include <Exceptions.hpp>
+#include <sstream>
+
+#include <rang.hpp>
+
+NormalError::NormalError(std::string_view error_msg)
+  : std::runtime_error([&] {
+      std::ostringstream os;
+      os << rang::style::bold << "Error:" << rang::style::reset << '\n' << error_msg << '\n';
+      return os.str();
+    }())
+{}
+
+UnexpectedError::UnexpectedError(std::string_view error_msg)
+  : std::runtime_error([&] {
+      std::ostringstream os;
+      os << rang::fgB::red << "Unexpected error:" << rang::style::reset << '\n' << error_msg << '\n';
+      return os.str();
+    }())
+{}
+
+NotImplementedError::NotImplementedError(std::string_view error_msg)
+  : std::runtime_error([&] {
+      std::ostringstream os;
+      os << rang::fgB::yellow << "Not implemented yet:" << rang::style::reset << '\n' << error_msg << '\n';
+      return os.str();
+    }())
+{}
diff --git a/src/utils/Exceptions.hpp b/src/utils/Exceptions.hpp
new file mode 100644
index 000000000..4256d02c8
--- /dev/null
+++ b/src/utils/Exceptions.hpp
@@ -0,0 +1,21 @@
+#ifndef EXCEPTIONS_H
+#define EXCEPTIONS_H
+
+#include <stdexcept>
+
+struct NormalError : public std::runtime_error
+{
+  NormalError(std::string_view error_msg);
+};
+
+struct UnexpectedError : public std::runtime_error
+{
+  UnexpectedError(std::string_view error_msg);
+};
+
+struct NotImplementedError : public std::runtime_error
+{
+  NotImplementedError(std::string_view error_msg);
+};
+
+#endif /* EXCEPTIONS_H */
diff --git a/src/utils/Messenger.cpp b/src/utils/Messenger.cpp
index 07eab1792..bb79a3719 100644
--- a/src/utils/Messenger.cpp
+++ b/src/utils/Messenger.cpp
@@ -12,8 +12,7 @@ Messenger::create(int& argc, char* argv[])
   if (Messenger::m_instance == nullptr) {
     Messenger::m_instance = new Messenger(argc, argv);
   } else {
-    std::cerr << "Messenger already created\n";
-    std::exit(1);
+    throw UnexpectedError("Messenger already created");
   }
 }
 
diff --git a/src/utils/Messenger.hpp b/src/utils/Messenger.hpp
index d3c3801be..85a95711a 100644
--- a/src/utils/Messenger.hpp
+++ b/src/utils/Messenger.hpp
@@ -18,6 +18,8 @@
 
 #include <PugsTraits.hpp>
 
+#include <Exceptions.hpp>
+
 #include <iostream>
 
 namespace parallel
@@ -232,8 +234,7 @@ class Messenger
       std::vector<MPI_Status> status_list(request_list.size());
       if (MPI_SUCCESS != MPI_Waitall(request_list.size(), &(request_list[0]), &(status_list[0]))) {
         // LCOV_EXCL_START
-        std::cerr << "Communication error!\n";
-        std::exit(1);
+        throw NormalError("Communication error");
         // LCOV_EXCL_STOP
       }
     }
diff --git a/src/utils/Partitioner.cpp b/src/utils/Partitioner.cpp
index 4bbee923c..2e22aa4f2 100644
--- a/src/utils/Partitioner.cpp
+++ b/src/utils/Partitioner.cpp
@@ -11,6 +11,8 @@
 #include <iostream>
 #include <vector>
 
+#include <Exceptions.hpp>
+
 Array<int>
 Partitioner::partition(const CSRGraph& graph)
 {
@@ -61,8 +63,7 @@ Partitioner::partition(const CSRGraph& graph)
       ParMETIS_V3_PartKway(&(vtxdist[0]), &(entries[0]), &(neighbors[0]), NULL, NULL, &wgtflag, &numflag, &ncon, &npart,
                            &(tpwgts[0]), &(ubvec[0]), &(options[0]), &edgecut, &(part[0]), &parmetis_comm);
     if (result == METIS_ERROR) {
-      std::cerr << "Metis Error\n";
-      std::exit(1);
+      throw UnexpectedError("Metis Error");
     }
 
     MPI_Comm_free(&parmetis_comm);
diff --git a/src/utils/PugsAssert.hpp b/src/utils/PugsAssert.hpp
index 1ca85e009..5c0a2cc62 100644
--- a/src/utils/PugsAssert.hpp
+++ b/src/utils/PugsAssert.hpp
@@ -7,13 +7,6 @@
 #include <rang.hpp>
 #include <string>
 
-template <typename ErrorType>
-void
-printAndThrow(const ErrorType& error)
-{
-  throw error;
-}
-
 class AssertError
 {
  private:
@@ -68,21 +61,21 @@ PRAGMA_DIAGNOSTIC_POP
 
 #else   // NDEBUG
 
-#define Assert(assertion, ...)                                                                       \
-  if (not _pugs_assert(assertion)) {                                                                 \
-    using vargs_t = decltype(std::make_tuple(__VA_ARGS__));                                          \
-    static_assert(std::tuple_size_v<vargs_t> <= 1, "too many arguments");                            \
-    if constexpr (std::tuple_size_v<vargs_t> == 0) {                                                 \
-      printAndThrow(AssertError(__FILE__, __LINE__, __PRETTY_FUNCTION__, #assertion));               \
-    } else {                                                                                         \
-      printAndThrow(AssertError(__FILE__, __LINE__, __PRETTY_FUNCTION__, #assertion, #__VA_ARGS__)); \
-    }                                                                                                \
+#define Assert(assertion, ...)                                                              \
+  if (not _pugs_assert(assertion)) {                                                        \
+    using vargs_t = decltype(std::make_tuple(__VA_ARGS__));                                 \
+    static_assert(std::tuple_size_v<vargs_t> <= 1, "too many arguments");                   \
+    if constexpr (std::tuple_size_v<vargs_t> == 0) {                                        \
+      throw AssertError(__FILE__, __LINE__, __PRETTY_FUNCTION__, #assertion);               \
+    } else {                                                                                \
+      throw AssertError(__FILE__, __LINE__, __PRETTY_FUNCTION__, #assertion, #__VA_ARGS__); \
+    }                                                                                       \
   }
 
 #endif   // NDEBUG
 
-// stores the current state of Assert macro. This is useful for instance to
-// noexcept management of Assert throws
+// store the current state of Assert macro. This is useful for
+// instance to noexcept management of Assert throws
 #ifdef NDEBUG
 #define NO_ASSERT true
 #else   // NDEBUG
diff --git a/src/utils/SignalManager.cpp b/src/utils/SignalManager.cpp
index c9cf877dc..a9981848d 100644
--- a/src/utils/SignalManager.cpp
+++ b/src/utils/SignalManager.cpp
@@ -81,6 +81,12 @@ SignalManager::handler(int signal)
       std::rethrow_exception(eptr);
     }
   }
+  catch (const NotImplementedError& not_implemented_error) {
+    std::cerr << not_implemented_error.what() << '\n';
+  }
+  catch (const UnexpectedError& unexpected_error) {
+    std::cerr << unexpected_error.what() << '\n';
+  }
   catch (const AssertError& assert_error) {
     std::cerr << assert_error << '\n';
   }
-- 
GitLab