Skip to content
Snippets Groups Projects
Commit 7882ebe0 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Add a mesh validity checker for Gmsh reader

It checks that
- no cell, face or edge is duplicated
- cells orientation is correct (positive volume)
parent 4f06d399
No related branches found
No related tags found
1 merge request!152Add a mesh validity checker for Gmsh reader
......@@ -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");
......
......@@ -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;
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment