MeshValidation.cpp
/**
* \file
* \author Karsten Rink
* \date 2013-04-04
* \brief Implementation of the MeshValidation class.
*
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#include "MeshValidation.h"
#include <numeric>
#include <algorithm>
#include <stack>
#include "BaseLib/Logging.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/Node.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/MeshEditing/MeshRevision.h"
#include "MeshLib/MeshSearch/NodeSearch.h"
#include "MeshLib/MeshSurfaceExtraction.h"
namespace MeshLib {
MeshValidation::MeshValidation(MeshLib::Mesh &mesh)
{
INFO ("Mesh Quality Control:");
INFO ("Looking for unused nodes...");
NodeSearch ns(mesh);
ns.searchUnused();
if (!ns.getSearchedNodeIDs().empty()) {
INFO("{:d} unused mesh nodes found.", ns.getSearchedNodeIDs().size());
}
MeshRevision rev(mesh);
INFO("Found {:d} potentially collapsible nodes.",
rev.getNumberOfCollapsableNodes());
const std::vector<ElementErrorCode> codes(
MeshLib::MeshValidation::testElementGeometry(mesh));
std::array<std::string,
static_cast<std::size_t>(ElementErrorFlag::MaxValue)>
output_str(MeshLib::MeshValidation::ElementErrorCodeOutput(codes));
for (auto& i : output_str)
{
INFO("{:s}", i);
}
}
std::vector<ElementErrorCode> MeshValidation::testElementGeometry(const MeshLib::Mesh &mesh, double min_volume)
{
INFO ("Testing mesh element geometry:");
const auto nErrorCodes(
static_cast<std::size_t>(ElementErrorFlag::MaxValue));
unsigned error_count[nErrorCodes];
std::fill_n(error_count, 4, 0);
const std::size_t nElements (mesh.getNumberOfElements());
const std::vector<MeshLib::Element*> &elements (mesh.getElements());
std::vector<ElementErrorCode> error_code_vector;
error_code_vector.reserve(nElements);
for (std::size_t i=0; i<nElements; ++i)
{
const ElementErrorCode e = elements[i]->validate();
error_code_vector.push_back(e);
if (e.none())
{
continue;
}
// increment error statistics
const std::bitset< static_cast<std::size_t>(ElementErrorFlag::MaxValue) > flags (static_cast< std::bitset<static_cast<std::size_t>(ElementErrorFlag::MaxValue)> >(e));
for (unsigned j = 0; j < nErrorCodes; ++j)
{
error_count[j] += flags[j];
}
}
// if a larger volume threshold is given, evaluate elements again to add them even if they are formally okay
if (min_volume > std::numeric_limits<double>::epsilon())
{
for (std::size_t i = 0; i < nElements; ++i)
{
if (elements[i]->getContent() < min_volume)
{
error_code_vector[i].set(ElementErrorFlag::ZeroVolume);
}
}
}
// output
const auto error_sum(static_cast<unsigned>(
std::accumulate(error_count, error_count + nErrorCodes, 0.0)));
if (error_sum != 0)
{
ElementErrorFlag flags[nErrorCodes] = { ElementErrorFlag::ZeroVolume, ElementErrorFlag::NonCoplanar,
ElementErrorFlag::NonConvex, ElementErrorFlag::NodeOrder };
for (std::size_t i = 0; i < nErrorCodes; ++i)
{
if (error_count[i])
{
INFO("{:d} elements found with {:s}.",
error_count[i],
ElementErrorCode::toString(flags[i]));
}
}
}
else
{
INFO ("No errors found.");
}
return error_code_vector;
}
std::array<std::string, static_cast<std::size_t>(ElementErrorFlag::MaxValue)>
MeshValidation::ElementErrorCodeOutput(const std::vector<ElementErrorCode> &error_codes)
{
const auto nErrorFlags(
static_cast<std::size_t>(ElementErrorFlag::MaxValue));
ElementErrorFlag flags[nErrorFlags] = {
ElementErrorFlag::ZeroVolume, ElementErrorFlag::NonCoplanar,
ElementErrorFlag::NonConvex, ElementErrorFlag::NodeOrder};
const std::size_t nElements(error_codes.size());
std::array<std::string,
static_cast<std::size_t>(ElementErrorFlag::MaxValue)>
output;
for (std::size_t i = 0; i < nErrorFlags; ++i)
{
unsigned count(0);
std::string elementIdStr;
for (std::size_t j = 0; j < nElements; ++j)
{
if (error_codes[j][flags[i]])
{
elementIdStr += (std::to_string(j) + ", ");
count++;
}
}
const std::string nErrorsStr = (count) ? std::to_string(count) : "No";
output[i] = (nErrorsStr + " elements found with " +
ElementErrorCode::toString(flags[i]) + ".\n");
if (count)
{
output[i] += ("ElementIDs: " + elementIdStr + "\n");
}
}
return output;
}
unsigned MeshValidation::detectHoles(MeshLib::Mesh const& mesh)
{
if (mesh.getDimension() == 1)
{
return 0;
}
auto boundary_mesh = MeshLib::BoundaryExtraction::getBoundaryElementsAsMesh(
mesh, "bulk_node_ids", "bulk_element_ids", "bulk_face_ids");
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);
auto 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());
}
// 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->getNumberOfNeighbors());
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