diff --git a/Applications/DataExplorer/mainwindow.cpp b/Applications/DataExplorer/mainwindow.cpp
index cce5464ee4f79c54fce2600c1e06c6596cbe62de..49b8307b8145b6d290c00462e8b66e16540c6bae 100644
--- a/Applications/DataExplorer/mainwindow.cpp
+++ b/Applications/DataExplorer/mainwindow.cpp
@@ -1145,76 +1145,90 @@ void MainWindow::callGMSH(std::vector<std::string>& selectedGeometries,
 
         if (!fileName.isEmpty())
         {
-            if (param4 == -1)
-            {  // adaptive meshing selected
-                FileIO::GMSH::GMSHInterface gmsh_io(
-                    _project.getGEOObjects(), true,
-                    FileIO::GMSH::MeshDensityAlgorithm::AdaptiveMeshDensity,
-                    param2, param3, param1, selectedGeometries, false, false);
-                BaseLib::IO::writeStringToFile(gmsh_io.writeToString(),
-                                               fileName.toStdString());
-            }
-            else
-            {  // homogeneous meshing selected
-                FileIO::GMSH::GMSHInterface gmsh_io(
-                    _project.getGEOObjects(), true,
-                    FileIO::GMSH::MeshDensityAlgorithm::FixedMeshDensity,
-                    param4, param3, param1, selectedGeometries, false, false);
-                BaseLib::IO::writeStringToFile(gmsh_io.writeToString(),
-                                               fileName.toStdString());
-            }
-
-            if (system(nullptr) != 0)  // command processor available
+            try
             {
-                QSettings settings;
-                std::string gmsh_path = settings.value("DataExplorerGmshPath")
-                                            .toString()
-                                            .toStdString();
+                if (param4 == -1)
+                {  // adaptive meshing selected
+                    FileIO::GMSH::GMSHInterface gmsh_io(
+                        _project.getGEOObjects(), true,
+                        FileIO::GMSH::MeshDensityAlgorithm::AdaptiveMeshDensity,
+                        param2, param3, param1, selectedGeometries, false,
+                        false);
+                    BaseLib::IO::writeStringToFile(gmsh_io.writeToString(),
+                                                   fileName.toStdString());
+                }
+                else
+                {  // homogeneous meshing selected
+                    FileIO::GMSH::GMSHInterface gmsh_io(
+                        _project.getGEOObjects(), true,
+                        FileIO::GMSH::MeshDensityAlgorithm::FixedMeshDensity,
+                        param4, param3, param1, selectedGeometries, false,
+                        false);
+                    BaseLib::IO::writeStringToFile(gmsh_io.writeToString(),
+                                                   fileName.toStdString());
+                }
 
-                if (!gmsh_path.empty())
+                if (system(nullptr) != 0)  // command processor available
                 {
-                    std::string fname(fileName.toStdString());
-                    std::string gmsh_command =
-                        "\"" + gmsh_path + "\" -2 -algo meshadapt " + fname;
-                    std::size_t pos(fname.rfind("."));
-                    if (pos != std::string::npos)
-                    {
-                        fname = fname.substr(0, pos);
-                    }
-                    gmsh_command += " -o " + fname + ".msh";
-                    // Newer gmsh versions write a newer file format for meshes
-                    // per default. At the moment we can't read this new format.
-                    // This is a switch for gmsh to write the 'old' file format.
-                    gmsh_command += " -format msh22";
-                    auto const return_value = std::system(gmsh_command.c_str());
-                    if (return_value != 0)
+                    QSettings settings;
+                    std::string gmsh_path =
+                        settings.value("DataExplorerGmshPath")
+                            .toString()
+                            .toStdString();
+
+                    if (!gmsh_path.empty())
                     {
-                        QString const message =
-                            "Execution of gmsh command returned non-zero "
-                            "status, " +
-                            QString::number(return_value);
-                        OGSError::box(message, "Error");
+                        std::string fname(fileName.toStdString());
+                        std::string gmsh_command =
+                            "\"" + gmsh_path + "\" -2 -algo meshadapt " + fname;
+                        std::size_t pos(fname.rfind("."));
+                        if (pos != std::string::npos)
+                        {
+                            fname = fname.substr(0, pos);
+                        }
+                        gmsh_command += " -o " + fname + ".msh";
+                        // Newer gmsh versions write a newer file format for
+                        // meshes per default. At the moment we can't read this
+                        // new format. This is a switch for gmsh to write the
+                        // 'old' file format.
+                        gmsh_command += " -format msh22";
+                        auto const return_value =
+                            std::system(gmsh_command.c_str());
+                        if (return_value != 0)
+                        {
+                            QString const message =
+                                "Execution of gmsh command returned non-zero "
+                                "status, " +
+                                QString::number(return_value);
+                            OGSError::box(message, "Error");
+                        }
+                        else
+                        {
+                            this->loadFile(ImportFileType::GMSH,
+                                           fileName.left(fileName.length() - 3)
+                                               .append("msh"));
+                        }
                     }
                     else
                     {
-                        this->loadFile(
-                            ImportFileType::GMSH,
-                            fileName.left(fileName.length() - 3).append("msh"));
+                        OGSError::box("Location of GMSH not specified.",
+                                      "Error");
                     }
                 }
                 else
                 {
-                    OGSError::box("Location of GMSH not specified.", "Error");
+                    OGSError::box(
+                        "Error executing command gmsh - no command processor "
+                        "available",
+                        "Error");
                 }
             }
-            else
+            catch (std::runtime_error& error)
             {
-                OGSError::box(
-                    "Error executing command gmsh - no command processor "
-                    "available",
-                    "Error");
+                OGSError::box(QString(error.what()) +
+                                  QString("\n Please cleanup the input data."),
+                              "ERROR");
             }
-
             if (delete_geo_file)
             {
                 BaseLib::removeFile(fileName.toStdString());
diff --git a/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp b/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
index 73b3236005d3f9bd236cc5bb308d3e09fbc2492c..9b72e148c8039da20e866cc1b29935763ca83d18 100644
--- a/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
+++ b/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
@@ -10,6 +10,8 @@
 
 #include "GMSHPolygonTree.h"
 
+#include <sstream>
+
 #include "GMSHAdaptiveMeshDensity.h"
 #include "GMSHFixedMeshDensity.h"
 #include "GeoLib/AnalyticalGeometry.h"
@@ -305,6 +307,8 @@ void GMSHPolygonTree::createGMSHPoints(std::vector<GMSHPoint*>& gmsh_pnts) const
     }
 
     const std::size_t n_plys(_plys.size());
+    std::stringstream error_messages;
+    error_messages.precision(std::numeric_limits<double>::digits10);
     for (std::size_t k(0); k < n_plys; k++)
     {
         const std::size_t n_pnts_in_ply(_plys[k]->getNumberOfPoints());
@@ -323,6 +327,22 @@ void GMSHPolygonTree::createGMSHPoints(std::vector<GMSHPoint*>& gmsh_pnts) const
                     *pnt, id,
                     _mesh_density_strategy.getMeshDensityAtPoint(pnt));
             }
+            else
+            {
+                auto const& p = *(_plys[k]->getPoint(j));
+                error_messages << "\n\tpoint with id " << p.getID()
+                               << " and coordinates (" << p[0] << ", " << p[1]
+                               << ", " << p[2]
+                               << ") is outside of the polygon.";
+            }
+        }
+    }
+    if (!parent())
+    {
+        auto const error_message = error_messages.str();
+        if (!error_message.empty())
+        {
+            OGS_FATAL(error_message);
         }
     }
 
diff --git a/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp b/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp
index 826aa39f357166ecdf2b960825e6586d91d2c066..5abbe995ebb11677f68daa171d4c08217f7d5e8e 100644
--- a/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp
+++ b/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp
@@ -57,6 +57,12 @@ int main(int argc, char* argv[])
     TCLAP::SwitchArg homogeneous_flag(
         "", "homogeneous", "Use Gmsh homogeneous meshing method.", false);
     cmd.add(homogeneous_flag);
+    TCLAP::ValueArg<std::string> merged_geometries_output(
+        "", "write_merged_geometries",
+        "file name for the output of the internal created geometry (*.gml) "
+        "(useful for debugging the data)",
+        false, "", "output file name");
+    cmd.add(merged_geometries_output);
 
     cmd.parse(argc, argv);
 
@@ -85,35 +91,57 @@ int main(int argc, char* argv[])
     bool const rotate = false;
     bool const keep_preprocessed_geometry = false;
 
-    if (homogeneous_flag.getValue())
-    {  // homogeneous meshing
-        double const average_mesh_density =
-            average_point_density_arg.getValue();
-        FileIO::GMSH::GMSHInterface gmsh_io(
-            geo_objects, true,
-            FileIO::GMSH::MeshDensityAlgorithm::FixedMeshDensity,
-            average_mesh_density, 0.0, 0, geo_names, rotate,
-            keep_preprocessed_geometry);
-        BaseLib::IO::writeStringToFile(gmsh_io.writeToString(),
-                                       geo_output_arg.getValue());
+    try
+    {
+        if (homogeneous_flag.getValue())
+        {  // homogeneous meshing
+            double const average_mesh_density =
+                average_point_density_arg.getValue();
+            FileIO::GMSH::GMSHInterface gmsh_io(
+                geo_objects, true,
+                FileIO::GMSH::MeshDensityAlgorithm::FixedMeshDensity,
+                average_mesh_density, 0.0, 0, geo_names, rotate,
+                keep_preprocessed_geometry);
+            BaseLib::IO::writeStringToFile(gmsh_io.writeToString(),
+                                           geo_output_arg.getValue());
+        }
+        else
+        {  // adaptive meshing
+            unsigned const max_number_of_points_in_quadtree_leaf =
+                max_number_of_points_in_quadtree_leaf_arg.getValue();
+            double const mesh_density_scaling_points =
+                mesh_density_scaling_points_arg.getValue();
+            double const mesh_density_scaling_stations =
+                mesh_density_scaling_stations_arg.getValue();
+            FileIO::GMSH::GMSHInterface gmsh_io(
+                geo_objects, true,
+                FileIO::GMSH::MeshDensityAlgorithm::AdaptiveMeshDensity,
+                mesh_density_scaling_points, mesh_density_scaling_stations,
+                max_number_of_points_in_quadtree_leaf, geo_names, rotate,
+                keep_preprocessed_geometry);
+            BaseLib::IO::writeStringToFile(gmsh_io.writeToString(),
+                                           geo_output_arg.getValue());
+        }
     }
-    else
-    {  // adaptive meshing
-        unsigned const max_number_of_points_in_quadtree_leaf =
-            max_number_of_points_in_quadtree_leaf_arg.getValue();
-        double const mesh_density_scaling_points =
-            mesh_density_scaling_points_arg.getValue();
-        double const mesh_density_scaling_stations =
-            mesh_density_scaling_stations_arg.getValue();
-        FileIO::GMSH::GMSHInterface gmsh_io(
-            geo_objects, true,
-            FileIO::GMSH::MeshDensityAlgorithm::AdaptiveMeshDensity,
-            mesh_density_scaling_points, mesh_density_scaling_stations,
-            max_number_of_points_in_quadtree_leaf, geo_names, rotate,
-            keep_preprocessed_geometry);
-        BaseLib::IO::writeStringToFile(gmsh_io.writeToString(),
-                                       geo_output_arg.getValue());
+    catch (std::runtime_error& error)
+    {
+        ERR("{}", error.what());
+        if (merged_geometries_output.isSet())
+        {
+            GeoLib::IO::BoostXmlGmlInterface xml(geo_objects);
+            xml.export_name = geo_objects.getGeometryNames().back();
+            BaseLib::IO::writeStringToFile(xml.writeToString(),
+                                           merged_geometries_output.getValue());
+        }
+        else
+        {
+            INFO(
+                "Hint: Using the command line flag "
+                "'--write_merged_geometries buggy_geometry.gml' allows "
+                "for better analysis of the problem.");
+        }
+        ERR("An error has occurred - programme will be terminated.");
+        return EXIT_FAILURE;
     }
-
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/Tests.cmake b/Applications/Utils/Tests.cmake
index 4d655385c69fc3bc324f492f838b92e918673002..4f8ff5827322422308367e1cdd1c4b6ef80ffeef 100644
--- a/Applications/Utils/Tests.cmake
+++ b/Applications/Utils/Tests.cmake
@@ -693,6 +693,23 @@ AddTest(
     square_1x1_homogeneous.geo
 )
 
+AddTest(
+    NAME LineIntersectingDomainBoundary
+    PATH MeshGeoToolsLib/geometryToGmshGeo/
+    WORKING_DIRECTORY ${Data_SOURCE_DIR}/MeshGeoToolsLib/geometryToGmshGeo
+    EXECUTABLE geometryToGmshGeo
+    EXECUTABLE_ARGS -i square_1x1.gml -i line_intersecting_square.gml -o ${Data_BINARY_DIR}/MeshGeoToolsLib/geometryToGmshGeo/square_1x1_with_intersecting_line.geo
+    REQUIREMENTS NOT OGS_USE_MPI
+)
+if (NOT OGS_USE_MPI)
+    set_tests_properties(
+        geometryToGmshGeo-LineIntersectingDomainBoundary
+        PROPERTIES
+            PASS_REGULAR_EXPRESSION
+            "ogs.*;ogs.*error.*\n\tpoint with id 5 and coordinates (1.001000000001, 0.6, 0) is outside of the polygon."
+    )
+endif()
+
 AddTest(
     NAME ResetPropertiesInPolygonalRegion_AllElementNodesInPolygon
     PATH MeshGeoToolsLib/ResetPropertiesInPolygonalRegion/
diff --git a/Tests/Data/MeshGeoToolsLib/geometryToGmshGeo/line_intersecting_square.gml b/Tests/Data/MeshGeoToolsLib/geometryToGmshGeo/line_intersecting_square.gml
new file mode 100644
index 0000000000000000000000000000000000000000..30182dee0ea3d41dd29bcf80424cea015d8bbd98
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/geometryToGmshGeo/line_intersecting_square.gml
@@ -0,0 +1,14 @@
+<?xml version="1.0" encoding="utf-8"?>
+<OpenGeoSysGLI>
+	<name>Line</name>
+	<points>
+		<point id="0" x="0.5" y="0.5" z="0"/>
+		<point id="1" x="1.001000000001" y="0.59999999999999998" z="0"/>
+	</points>
+	<polylines>
+		<polyline id="0" name="Line">
+			<pnt>0</pnt>
+			<pnt>1</pnt>
+		</polyline>
+	</polylines>
+</OpenGeoSysGLI>