diff --git a/Applications/DataExplorer/DataView/MeshAnalysis.ui b/Applications/DataExplorer/DataView/MeshAnalysis.ui index 70db2547006ace3442b66b931a3deff2d9c7c12a..d97b90e66a9126b1ab2a261218f0f336b1f98320 100644 --- a/Applications/DataExplorer/DataView/MeshAnalysis.ui +++ b/Applications/DataExplorer/DataView/MeshAnalysis.ui @@ -7,7 +7,7 @@ <x>0</x> <y>0</y> <width>483</width> - <height>650</height> + <height>787</height> </rect> </property> <property name="windowTitle"> @@ -348,6 +348,13 @@ </property> </widget> </item> + <item row="5" column="0" colspan="4"> + <widget class="QLabel" name="meshHoleOutputLabel"> + <property name="text"> + <string/> + </property> + </widget> + </item> </layout> </widget> </item> diff --git a/Applications/DataExplorer/DataView/MeshAnalysisDialog.cpp b/Applications/DataExplorer/DataView/MeshAnalysisDialog.cpp index a3986619a39eb7bdc5980df6328c7542a204f804..05a63c4df3a9928f036b09cc0c3d7a31238bd0c8 100644 --- a/Applications/DataExplorer/DataView/MeshAnalysisDialog.cpp +++ b/Applications/DataExplorer/DataView/MeshAnalysisDialog.cpp @@ -60,6 +60,10 @@ void MeshAnalysisDialog::on_startButton_pressed() *mesh, this->zeroVolumeThreshold->text().toDouble() + std::numeric_limits<double>::epsilon())); this->elementsGroupBox->setTitle("Elements (out of " + QString::number(mesh->getNElements()) + ")"); this->elementsMsgOutput(element_error_codes); + + unsigned const n_holes (MeshLib::MeshValidation::detectHoles(*mesh)); + if (n_holes>0) + this->meshHoleOutputLabel->setText("<strong>" + QString::number(n_holes) + " hole(s) found within the mesh</strong>"); } void MeshAnalysisDialog::nodesMsgOutput(std::vector<std::size_t> const& node_ids, std::vector<std::size_t> const& collapsibleNodeIds) diff --git a/Applications/Utils/MeshEdit/checkMesh.cpp b/Applications/Utils/MeshEdit/checkMesh.cpp index 5fe055778bcbdded6c5cdd3a2bbd95eefa4e43d6..0314bf8647b48c40a22f313563e123e2295605f0 100644 --- a/Applications/Utils/MeshEdit/checkMesh.cpp +++ b/Applications/Utils/MeshEdit/checkMesh.cpp @@ -91,8 +91,16 @@ int main(int argc, char *argv[]) // Validation MeshLib::MeshValidation validation(*const_cast<MeshLib::Mesh*>(mesh)); // MeshValidation outputs error messages /* Remark: MeshValidation can modify the original mesh */ - } + unsigned const n_holes (MeshLib::MeshValidation::detectHoles(*mesh)); + if (n_holes > 0) + { + INFO("%d hole(s) detected within the mesh", n_holes); + } + else + INFO ("No holes found within the mesh."); + } + delete mesh; delete custom_format; delete logog_cout; diff --git a/MeshLib/MeshQuality/MeshValidation.cpp b/MeshLib/MeshQuality/MeshValidation.cpp index 15dc4327a0d202f70a18a55d43fbe92ad71fad99..940461ae573d9995a9daeaca09dbf2fcf4d36912 100644 --- a/MeshLib/MeshQuality/MeshValidation.cpp +++ b/MeshLib/MeshQuality/MeshValidation.cpp @@ -15,6 +15,8 @@ #include "MeshValidation.h" #include <numeric> +#include <algorithm> +#include <stack> #include "logog/include/logog.hpp" @@ -23,6 +25,7 @@ #include "MeshLib/Elements/Element.h" #include "MeshLib/MeshEditing/removeMeshNodes.h" #include "MeshLib/MeshEditing/MeshRevision.h" +#include "MeshLib/MeshSurfaceExtraction.h" namespace MeshLib { @@ -147,4 +150,44 @@ MeshValidation::ElementErrorCodeOutput(const std::vector<ElementErrorCode> &erro return output; } +unsigned MeshValidation::detectHoles(MeshLib::Mesh const& mesh) +{ + MeshLib::Mesh* boundary_mesh (MeshSurfaceExtraction::getMeshBoundary(mesh)); + std::vector<MeshLib::Element*> const& elements (boundary_mesh->getElements()); + + std::vector<unsigned> sfc_idx (elements.size(), std::numeric_limits<unsigned>::max()); + unsigned current_surface_id (0); + std::vector<unsigned>::const_iterator it = sfc_idx.cbegin(); + + while (it != sfc_idx.cend()) + { + std::size_t const idx = static_cast<std::size_t>(std::distance(sfc_idx.cbegin(), it)); + trackSurface(elements[idx], sfc_idx, current_surface_id++); + it = std::find(sfc_idx.cbegin(), sfc_idx.cend(), std::numeric_limits<unsigned>::max()); + } + delete boundary_mesh; + + // Subtract "1" from the number of surfaces found to get the number of holes. + return (--current_surface_id); +} + +void MeshValidation::trackSurface(MeshLib::Element const* element, std::vector<unsigned> &sfc_idx, unsigned const current_index) +{ + std::stack<MeshLib::Element const*> elem_stack; + elem_stack.push(element); + while (!elem_stack.empty()) + { + MeshLib::Element const*const elem = elem_stack.top(); + elem_stack.pop(); + sfc_idx[elem->getID()] = current_index; + std::size_t const n_neighbors (elem->getNNeighbors()); + for (std::size_t i=0; i<n_neighbors; ++i) + { + MeshLib::Element const* neighbor (elem->getNeighbor(i)); + if ( neighbor != nullptr && sfc_idx[neighbor->getID()] == std::numeric_limits<unsigned>::max()) + elem_stack.push(neighbor); + } + } +} + } // end namespace MeshLib diff --git a/MeshLib/MeshQuality/MeshValidation.h b/MeshLib/MeshQuality/MeshValidation.h index 4d9905e447bb26a9f79b828c22603bcf2d2565fe..5dc222a0b1e5c1dcf466b2e41edfda14100c783f 100644 --- a/MeshLib/MeshQuality/MeshValidation.h +++ b/MeshLib/MeshQuality/MeshValidation.h @@ -23,6 +23,7 @@ namespace MeshLib { class Mesh; + class Element; /** * \brief A collection of methods for testing mesh quality and correctness @@ -65,7 +66,26 @@ public: static std::array<std::string, static_cast<std::size_t>(ElementErrorFlag::MaxValue)> ElementErrorCodeOutput(const std::vector<ElementErrorCode> &error_codes); + /** + * Tests if holes are located within the mesh. + * In this context, a hole is a boundary of an element with no neighbor that cannot be reached from + * the actual boundary of the mesh. Examples include a missing triangle in a 2D mesh or a missing + * Tetrahedron in a 3D mesh. + * Note, that this method does not work when complex 3D structures are build from 2D mesh elements, + * e.g. using the LayeredVolume-class, where more than two 2D elements may share an edge. + * @param mesh The mesh that is tested + * @return The number of holes that have been found. + */ + static unsigned detectHoles(MeshLib::Mesh const& mesh); + private: + /** Finds all surface elements that can be reached from element. All elements that are found in this + * way are marked in the global sfc_idx vector using the current_index. + * @param element The mesh element from which the search is started + * @param sfc_idx The global index vector notifying to which surface elements belong + * @param current_index The index that all elements reachable from element will be assigned in sfc_idx + */ + static void trackSurface(MeshLib::Element const*const element, std::vector<unsigned> &sfc_idx, unsigned const current_index); }; diff --git a/MeshLib/MeshSurfaceExtraction.cpp b/MeshLib/MeshSurfaceExtraction.cpp index 6b81367bc0e8c95d24b1405e24cbfe944fceabfe..1b3e8c58cff59a8e0bb9df68af55931f9f8a0042 100644 --- a/MeshLib/MeshSurfaceExtraction.cpp +++ b/MeshLib/MeshSurfaceExtraction.cpp @@ -20,11 +20,14 @@ #include "GeoLib/PointWithID.h" -#include "Mesh.h" -#include "Elements/Face.h" -#include "Elements/Cell.h" -#include "Elements/Tri.h" -#include "Elements/Quad.h" +#include "MeshLib/Mesh.h" +#include "MeshLib/Elements/Line.h" +#include "MeshLib/Elements/Face.h" +#include "MeshLib/Elements/Cell.h" +#include "MeshLib/Elements/Tri.h" +#include "MeshLib/Elements/Quad.h" +#include "MeshLib/MeshEditing/DuplicateMeshComponents.h" +#include "MeshLib/MeshQuality/MeshValidation.h" namespace MeshLib { @@ -118,6 +121,40 @@ MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh, return result; } +MeshLib::Mesh* MeshSurfaceExtraction::getMeshBoundary(const MeshLib::Mesh &mesh) +{ + if (mesh.getDimension()==1) + return nullptr; + + // For 3D meshes return the 2D surface + if (mesh.getDimension()==3) + { + MathLib::Vector3 dir(0,0,0); + return getMeshSurface(mesh, dir, 90); + } + + // For 2D meshes return the boundary lines + std::vector<MeshLib::Node*> nodes = MeshLib::copyNodeVector(mesh.getNodes()); + std::vector<MeshLib::Element*> boundary_elements; + + std::vector<MeshLib::Element*> const& org_elems (mesh.getElements()); + for (auto it=org_elems.begin(); it!=org_elems.end(); ++it) + { + MeshLib::Element* elem (*it); + std::size_t const n_edges (elem->getNEdges()); + for (std::size_t i=0; i<n_edges; ++i) + if (elem->getNeighbor(i) == nullptr) + { + MeshLib::Element const*const edge (elem->getEdge(i)); + boundary_elements.push_back(MeshLib::copyElement(edge, nodes)); + delete edge; + } + } + MeshLib::Mesh* result = new MeshLib::Mesh("Boundary Mesh", nodes, boundary_elements); + MeshLib::MeshValidation::removeUnusedMeshNodes(*result); + return result; +} + void MeshSurfaceExtraction::get2DSurfaceElements(const std::vector<MeshLib::Element*> &all_elements, std::vector<MeshLib::Element*> &sfc_elements, const MathLib::Vector3 &dir, double angle, unsigned mesh_dimension) { if (mesh_dimension<2 || mesh_dimension>3) diff --git a/MeshLib/MeshSurfaceExtraction.h b/MeshLib/MeshSurfaceExtraction.h index 73d7c4cb046b507741fb27b43271378e65c6567a..8e634f02de5e120613c22e81f5b09cbb2bcd77b7 100644 --- a/MeshLib/MeshSurfaceExtraction.h +++ b/MeshLib/MeshSurfaceExtraction.h @@ -39,11 +39,11 @@ public: /// Returns a vector of the areas assigned to each node on a surface mesh. static std::vector<double> getSurfaceAreaForNodes(const MeshLib::Mesh &mesh); - /// Returns the surface nodes of a layered mesh. + /// Returns the surface nodes of a mesh. static std::vector<GeoLib::PointWithID*> getSurfaceNodes(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, double angle); /** - * Returns the 2d-element mesh representing the surface of the given layered mesh. + * Returns the 2d-element mesh representing the surface of the given mesh. * \param mesh The original mesh * \param dir The direction in which face normals have to point to be considered surface elements * \param angle The angle of the allowed deviation from the given direction (0 <= angle <= 90 degrees) @@ -52,6 +52,15 @@ public: */ static MeshLib::Mesh* getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, double angle, bool keepOriginalNodeIds = false); + /** + * Returns the boundary of mesh, i.e. lines for 2D meshes and surfaces for 3D meshes. + * Note, that this method also returns inner boundaries and might give unexpected results + * when the mesh geometry is not strict (e.g. more than two triangles sharing an edge). + * \param mesh The original mesh of dimension d + * \return A mesh of dimension (d-1) representing the boundary of the mesh. + */ + static MeshLib::Mesh* getMeshBoundary(const MeshLib::Mesh &mesh); + private: /// Functionality needed for getSurfaceNodes() and getMeshSurface() static void get2DSurfaceElements(const std::vector<MeshLib::Element*> &all_elements, std::vector<MeshLib::Element*> &sfc_elements, const MathLib::Vector3 &dir, double angle, unsigned mesh_dimension); diff --git a/Tests/MeshLib/TestMeshValidation.cpp b/Tests/MeshLib/TestMeshValidation.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ecc02ee6cbb8ac9c00345a489bd1800514986e41 --- /dev/null +++ b/Tests/MeshLib/TestMeshValidation.cpp @@ -0,0 +1,128 @@ +/** + * @file TestMeshValidation.cpp + * @author Karsten Rink + * @date 2015-01-29 + * @brief Tests for MeshRevision class + * + * @copyright + * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/LICENSE.txt + */ + +#include "gtest/gtest.h" + +#include "MeshLib/Mesh.h" +#include "MeshLib/Node.h" +#include "MeshLib/MeshQuality/MeshValidation.h" +#include "MeshLib/MeshGenerators/MeshGenerator.h" +#include "MeshLib/MeshGenerators/ConvertRasterToMesh.h" +#include "MeshLib/MeshEditing/DuplicateMeshComponents.h" +#include "MeshLib/Elements/Element.h" + +#include "GeoLib/Raster.h" + +TEST(MeshValidation, UnusedNodes) +{ + std::array<double, 12> pix = {0,0.1,0.2,0.1,0,0,0.1,0,0,0,-0.1,0}; + GeoLib::Raster raster(4,3,0,0,1,pix.begin(), pix.end()); + MeshLib::ConvertRasterToMesh conv(raster, MeshElemType::TRIANGLE, MeshLib::UseIntensityAs::ELEVATION); + MeshLib::Mesh* mesh = conv.execute(); + std::vector<std::size_t> u_nodes = MeshLib::MeshValidation::findUnusedMeshNodes(*mesh); + ASSERT_EQ(0, u_nodes.size()); + + std::vector<MeshLib::Node*> nodes = MeshLib::copyNodeVector(mesh->getNodes()); + nodes.push_back(new MeshLib::Node(-1,-1,-1)); + std::vector<MeshLib::Element*> elems = MeshLib::copyElementVector(mesh->getElements(),nodes); + MeshLib::Mesh mesh2("mesh2", nodes, elems); + u_nodes = MeshLib::MeshValidation::findUnusedMeshNodes(mesh2); + ASSERT_EQ(1, u_nodes.size()); + ASSERT_EQ(nodes.back()->getID(), u_nodes[0]); + + delete mesh; +} + +TEST(MeshValidation, DetectHolesTri) +{ + std::array<double, 12> pix = {0,0.1,0.2,0.1,0,0,0.1,0,0,0,-0.1,0}; + GeoLib::Raster raster(4,3,0,0,1,pix.begin(), pix.end()); + MeshLib::ConvertRasterToMesh conv(raster, MeshElemType::TRIANGLE, MeshLib::UseIntensityAs::ELEVATION); + MeshLib::Mesh* mesh = conv.execute(); + unsigned n_holes = MeshLib::MeshValidation::detectHoles(*mesh); + ASSERT_EQ(0, n_holes); + + { + std::vector<MeshLib::Node*> nodes = MeshLib::copyNodeVector(mesh->getNodes()); + std::vector<MeshLib::Element*> elems = MeshLib::copyElementVector(mesh->getElements(),nodes); + elems.erase(elems.begin()+12); + MeshLib::Mesh mesh2("mesh2", nodes, elems); + n_holes = MeshLib::MeshValidation::detectHoles(mesh2); + ASSERT_EQ(1, n_holes); + } + + { + std::vector<MeshLib::Node*> nodes = MeshLib::copyNodeVector(mesh->getNodes()); + std::vector<MeshLib::Element*> elems = MeshLib::copyElementVector(mesh->getElements(),nodes); + elems.erase(elems.begin()+11); + elems.erase(elems.begin()+11); + MeshLib::Mesh mesh2("mesh2", nodes, elems); + n_holes = MeshLib::MeshValidation::detectHoles(mesh2); + ASSERT_EQ(1, n_holes); + } + + { + std::vector<MeshLib::Node*> nodes = MeshLib::copyNodeVector(mesh->getNodes()); + std::vector<MeshLib::Element*> elems = MeshLib::copyElementVector(mesh->getElements(),nodes); + elems.erase(elems.begin()+10); + elems.erase(elems.begin()+12); + MeshLib::Mesh mesh2("mesh2", nodes, elems); + n_holes = MeshLib::MeshValidation::detectHoles(mesh2); + ASSERT_EQ(2, n_holes); + } + + delete mesh; +} + +TEST(MeshValidation, DetectHolesHex) +{ + std::size_t xs=5, ys=4, zs=4; + GeoLib::Point origin(0,0,0); + MeshLib::Mesh* mesh = MeshLib::MeshGenerator::generateRegularHexMesh(5,4,4,1,1,1,origin, "mesh"); + unsigned n_holes = MeshLib::MeshValidation::detectHoles(*mesh); + ASSERT_EQ(0, n_holes); + + { + std::vector<MeshLib::Node*> nodes = MeshLib::copyNodeVector(mesh->getNodes()); + std::vector<MeshLib::Element*> elems = MeshLib::copyElementVector(mesh->getElements(),nodes); + elems.erase(elems.begin()+27); + MeshLib::Mesh mesh2("mesh2", nodes, elems); + n_holes = MeshLib::MeshValidation::detectHoles(mesh2); + ASSERT_EQ(1, n_holes); + } + + { + std::vector<MeshLib::Node*> nodes = MeshLib::copyNodeVector(mesh->getNodes()); + std::vector<MeshLib::Element*> elems = MeshLib::copyElementVector(mesh->getElements(),nodes); + elems.erase(elems.begin()+28); + elems.erase(elems.begin()+27); + MeshLib::Mesh mesh2("mesh2", nodes, elems); + n_holes = MeshLib::MeshValidation::detectHoles(mesh2); + ASSERT_EQ(1, n_holes); + } + + { + std::vector<MeshLib::Node*> nodes = MeshLib::copyNodeVector(mesh->getNodes()); + std::vector<MeshLib::Element*> elems = MeshLib::copyElementVector(mesh->getElements(),nodes); + elems.erase(elems.begin()+29); + elems.erase(elems.begin()+27); + MeshLib::Mesh mesh2("mesh2", nodes, elems); + n_holes = MeshLib::MeshValidation::detectHoles(mesh2); + ASSERT_EQ(1, n_holes); + } + delete mesh; +} + + + +