diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 2420eb4dba02662207ec8e0e15e2300c36d72ab2..e555ca224c1d02f11e11696542774f2e470b6432 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -1516,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);
@@ -1535,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);
@@ -1554,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");
diff --git a/src/mesh/MeshBuilderBase.cpp b/src/mesh/MeshBuilderBase.cpp
index 9e62baf56e7915379d65a5c88c21fe7411fca253..a3a8440bd25499d8aea34a08e7c836b66a4c3a32 100644
--- a/src/mesh/MeshBuilderBase.cpp
+++ b/src/mesh/MeshBuilderBase.cpp
@@ -5,12 +5,14 @@
#include <mesh/ConnectivityDispatcher.hpp>
#include <mesh/ItemId.hpp>
#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
#include <vector>
-template <int Dimension>
+template <size_t Dimension>
void
MeshBuilderBase::_dispatch()
{
@@ -28,8 +30,8 @@ MeshBuilderBase::_dispatch()
NodeValue<Rd> xr;
m_mesh = std::make_shared<MeshType>(connectivity, xr);
}
- const MeshType& mesh = static_cast<const MeshType&>(*m_mesh);
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*m_mesh);
ConnectivityDispatcher<Dimension> dispatcher(mesh.connectivity());
std::shared_ptr dispatched_connectivity = dispatcher.dispatchedConnectivity();
@@ -41,3 +43,165 @@ MeshBuilderBase::_dispatch()
template void MeshBuilderBase::_dispatch<1>();
template void MeshBuilderBase::_dispatch<2>();
template void MeshBuilderBase::_dispatch<3>();
+
+template <size_t Dimension>
+void
+MeshBuilderBase::_checkMesh() const
+{
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+
+ if (not m_mesh) {
+ throw UnexpectedError("mesh is not built yet");
+ }
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*m_mesh);
+
+ const ConnectivityType& connectivity = mesh.connectivity();
+
+ if constexpr (Dimension > 2) { // check for duplicated edges
+ auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
+
+ std::vector<std::vector<EdgeId>> node_edges(mesh.numberOfNodes());
+ for (EdgeId edge_id = 0; edge_id < mesh.numberOfEdges(); ++edge_id) {
+ auto node_list = edge_to_node_matrix[edge_id];
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ node_edges[node_list[i_node]].push_back(edge_id);
+ }
+ }
+
+ for (auto&& edge_list : node_edges) {
+ std::sort(edge_list.begin(), edge_list.end());
+ }
+
+ for (EdgeId edge_id = 0; edge_id < mesh.numberOfEdges(); ++edge_id) {
+ auto node_list = edge_to_node_matrix[edge_id];
+ std::vector<EdgeId> intersection = node_edges[node_list[0]];
+ for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
+ std::vector<EdgeId> local_intersection;
+ std::set_intersection(intersection.begin(), intersection.end(), //
+ node_edges[node_list[i_node]].begin(), node_edges[node_list[i_node]].end(),
+ std::back_inserter(local_intersection));
+ std::swap(local_intersection, intersection);
+ if (intersection.size() < 2) {
+ break;
+ }
+ }
+
+ if (intersection.size() > 1) {
+ std::ostringstream error_msg;
+ error_msg << "invalid mesh.\n\tFollowing edges\n";
+ for (EdgeId edge_id : intersection) {
+ error_msg << "\t - id=" << edge_id << " number=" << connectivity.edgeNumber()[edge_id] << '\n';
+ }
+ error_msg << "\tare duplicated";
+ throw NormalError(error_msg.str());
+ }
+ }
+ }
+
+ if constexpr (Dimension > 1) { // check for duplicated faces
+ auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+
+ std::vector<std::vector<FaceId>> node_faces(mesh.numberOfNodes());
+ for (FaceId face_id = 0; face_id < mesh.numberOfFaces(); ++face_id) {
+ auto node_list = face_to_node_matrix[face_id];
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ node_faces[node_list[i_node]].push_back(face_id);
+ }
+ }
+
+ for (auto&& face_list : node_faces) {
+ std::sort(face_list.begin(), face_list.end());
+ }
+
+ for (FaceId face_id = 0; face_id < mesh.numberOfFaces(); ++face_id) {
+ auto node_list = face_to_node_matrix[face_id];
+ std::vector<FaceId> intersection = node_faces[node_list[0]];
+ for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
+ std::vector<FaceId> local_intersection;
+ std::set_intersection(intersection.begin(), intersection.end(), //
+ node_faces[node_list[i_node]].begin(), node_faces[node_list[i_node]].end(),
+ std::back_inserter(local_intersection));
+ std::swap(local_intersection, intersection);
+ if (intersection.size() < 2) {
+ break;
+ }
+ }
+
+ if (intersection.size() > 1) {
+ std::ostringstream error_msg;
+ error_msg << "invalid mesh.\n\tFollowing faces\n";
+ for (FaceId face_id : intersection) {
+ error_msg << "\t - id=" << face_id << " number=" << connectivity.faceNumber()[face_id] << '\n';
+ }
+ error_msg << "\tare duplicated";
+ throw NormalError(error_msg.str());
+ }
+ }
+ }
+
+ auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+
+ { // check for duplicated cells
+ std::vector<std::vector<CellId>> node_cells(mesh.numberOfNodes());
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto node_list = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ node_cells[node_list[i_node]].push_back(cell_id);
+ }
+ }
+
+ for (auto&& cell_list : node_cells) {
+ std::sort(cell_list.begin(), cell_list.end());
+ }
+
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto node_list = cell_to_node_matrix[cell_id];
+ std::vector<CellId> intersection = node_cells[node_list[0]];
+ for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
+ std::vector<CellId> local_intersection;
+ std::set_intersection(intersection.begin(), intersection.end(), //
+ node_cells[node_list[i_node]].begin(), node_cells[node_list[i_node]].end(),
+ std::back_inserter(local_intersection));
+ std::swap(local_intersection, intersection);
+ if (intersection.size() < 2) {
+ break;
+ }
+ }
+
+ if (intersection.size() > 1) {
+ std::ostringstream error_msg;
+ error_msg << "invalid mesh.\n\tFollowing cells\n";
+ for (CellId cell_id : intersection) {
+ error_msg << "\t - id=" << cell_id << " number=" << connectivity.cellNumber()[cell_id] << '\n';
+ }
+ error_msg << "\tare duplicated";
+ throw NormalError(error_msg.str());
+ }
+ }
+ }
+
+ const auto& Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+ const auto& xr = mesh.xr();
+
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ double cell_volume = 0;
+ auto cell_nodes = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ cell_volume += dot(Cjr(cell_id, i_node), xr[cell_nodes[i_node]]);
+ }
+
+ if (cell_volume <= 0) {
+ std::ostringstream error_msg;
+ error_msg << "invalid mesh.\n\tThe following cell\n";
+ error_msg << "\t - id=" << cell_id << " number=" << connectivity.cellNumber()[cell_id] << '\n';
+ error_msg << "\thas non-positive volume: " << cell_volume / Dimension;
+ throw NormalError(error_msg.str());
+ }
+ }
+}
+
+template void MeshBuilderBase::_checkMesh<1>() const;
+template void MeshBuilderBase::_checkMesh<2>() const;
+template void MeshBuilderBase::_checkMesh<3>() const;
diff --git a/src/mesh/MeshBuilderBase.hpp b/src/mesh/MeshBuilderBase.hpp
index 6c0328bb8ee624e19d42fb4660ba3665ff7164de..63c55ec565b1f48023906f5aa602db13ef76d62c 100644
--- a/src/mesh/MeshBuilderBase.hpp
+++ b/src/mesh/MeshBuilderBase.hpp
@@ -10,9 +10,12 @@ class MeshBuilderBase
protected:
std::shared_ptr<const IMesh> m_mesh;
- template <int Dimension>
+ template <size_t Dimension>
void _dispatch();
+ template <size_t Dimension>
+ void _checkMesh() const;
+
public:
std::shared_ptr<const IMesh>
mesh() const