Add a mesh validity checker for Gmsh reader

It checks that

• no cell, face or edge is duplicated
• cells orientation is correct (positive volume)

src/mesh/GmshReader.cpp

@@ -70,7 +70,6 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_point_list.second.size(); ++j) {
       point_list[j] = ref_point_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_point_list.first);
 
     bool is_boundary = true;
     for (size_t i_node = 0; i_node < point_list.size(); ++i_node) {
@@ -80,7 +79,16 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
       }
     }
 
-    descriptor.addRefItemList(RefNodeList(physical_ref_id.refId(), point_list, is_boundary));
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_point_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_point_list.first};
+      }
+    }();
+
+    descriptor.addRefItemList(RefNodeList(ref_id, point_list, is_boundary));
   }
 
   std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
@@ -95,8 +103,17 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_cell_list.second.size(); ++j) {
       cell_list[j] = ref_cell_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_cell_list.first);
-    descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list, false));
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_cell_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_cell_list.first};
+      }
+    }();
+
+    descriptor.addRefItemList(RefCellList(ref_id, cell_list, false));
   }
 
   descriptor.cell_owner_vector.resize(nb_cells);
@@ -157,8 +174,17 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_cell_list.second.size(); ++j) {
       cell_list[j] = ref_cell_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_cell_list.first);
-    descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list, false));
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_cell_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_cell_list.first};
+      }
+    }();
+
+    descriptor.addRefItemList(RefCellList(ref_id, cell_list, false));
   }
 
   ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
@@ -229,7 +255,15 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_face_list.second.size(); ++j) {
       face_list[j] = ref_face_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_face_list.first);
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_face_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_face_list.first};
+      }
+    }();
 
     bool is_boundary = true;
     for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
@@ -239,7 +273,7 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
       }
     }
 
-    descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list, is_boundary});
+    descriptor.addRefItemList(RefFaceList{ref_id, face_list, is_boundary});
   }
 
   Array<bool> is_boundary_node(descriptor.node_number_vector.size());
@@ -265,7 +299,15 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_point_list.second.size(); ++j) {
       point_list[j] = ref_point_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_point_list.first);
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_point_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_point_list.first};
+      }
+    }();
 
     bool is_boundary = true;
     for (size_t i_node = 0; i_node < point_list.size(); ++i_node) {
@@ -274,7 +316,7 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
       }
     }
 
-    descriptor.addRefItemList(RefNodeList(physical_ref_id.refId(), point_list, is_boundary));
+    descriptor.addRefItemList(RefNodeList(ref_id, point_list, is_boundary));
   }
 
   descriptor.cell_owner_vector.resize(nb_cells);
@@ -372,8 +414,17 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_cell_list.second.size(); ++j) {
       cell_list[j] = ref_cell_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_cell_list.first);
-    descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list, false));
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_cell_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_cell_list.first};
+      }
+    }();
+
+    descriptor.addRefItemList(RefCellList(ref_id, cell_list, false));
   }
 
   ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
@@ -481,7 +532,15 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
       for (size_t j = 0; j < ref_face_list.second.size(); ++j) {
         face_list[j] = ref_face_list.second[j];
       }
-      const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_face_list.first);
+
+      auto ref_id = [&] {
+        if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_face_list.first);
+            i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+          return i_physical_ref->second.refId();
+        } else {
+          return RefId{ref_face_list.first};
+        }
+      }();
 
       bool is_boundary = true;
       for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
@@ -490,8 +549,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
           break;
         }
       }
-
-      descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list, is_boundary});
+      descriptor.addRefItemList(RefFaceList(ref_id, face_list, is_boundary));
     }
   }
 
@@ -566,7 +624,15 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
       for (size_t j = 0; j < ref_edge_list.second.size(); ++j) {
         edge_list[j] = ref_edge_list.second[j];
       }
-      const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_edge_list.first);
+
+      auto ref_id = [&] {
+        if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_edge_list.first);
+            i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+          return i_physical_ref->second.refId();
+        } else {
+          return RefId{ref_edge_list.first};
+        }
+      }();
 
       bool is_boundary = true;
       for (size_t i_node = 0; i_node < edge_list.size(); ++i_node) {
@@ -575,7 +641,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
         }
       }
 
-      descriptor.addRefItemList(RefEdgeList{physical_ref_id.refId(), edge_list, is_boundary});
+      descriptor.addRefItemList(RefEdgeList(ref_id, edge_list, is_boundary));
     }
   }
 
@@ -602,7 +668,15 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_point_list.second.size(); ++j) {
       point_list[j] = ref_point_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_point_list.first);
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_point_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_point_list.first};
+      }
+    }();
 
     bool is_boundary = true;
     for (size_t i_node = 0; i_node < point_list.size(); ++i_node) {
@@ -611,7 +685,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
       }
     }
 
-    descriptor.addRefItemList(RefNodeList(physical_ref_id.refId(), point_list, is_boundary));
+    descriptor.addRefItemList(RefNodeList(ref_id, point_list, is_boundary));
   }
 
   descriptor.cell_owner_vector.resize(nb_cells);
@@ -740,6 +814,7 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
     this->__proceedData();
   }
   std::cout << std::flush;
+
   if (parallel::size() > 1) {
     std::cout << "Sequential mesh read! Need to be dispatched\n" << std::flush;
 
@@ -763,6 +838,7 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
 
       return *begin(dimension_set);
     }();
+
     switch (mesh_dimension) {
     case 1: {
       this->_dispatch<1>();
@@ -1440,6 +1516,7 @@ GmshReader::__proceedData()
       xr[i][2] = m_mesh_data.__vertices[i][2];
     }
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
+    this->_checkMesh<3>();
 
   } else if (dot(dimension2_mask, elementNumber) > 0) {
     const size_t nb_cells = dot(dimension2_mask, elementNumber);
@@ -1459,6 +1536,7 @@ GmshReader::__proceedData()
       xr[i][1] = m_mesh_data.__vertices[i][1];
     }
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
+    this->_checkMesh<2>();
 
   } else if (dot(dimension1_mask, elementNumber) > 0) {
     const size_t nb_cells = dot(dimension1_mask, elementNumber);
@@ -1478,6 +1556,7 @@ GmshReader::__proceedData()
     }
 
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
+    this->_checkMesh<1>();
 
   } else {
     throw NormalError("using a dimension 0 mesh is forbidden");