MeshValidation.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include "MeshValidation.h"
 
#include <algorithm>
#include <numeric>
#include <stack>
 
#include "BaseLib/Logging.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshSearch/NodeSearch.h"
#include "MeshLib/Node.h"
#include "MeshToolsLib/MeshEditing/MeshRevision.h"
#include "MeshToolsLib/MeshSurfaceExtraction.h"
 
namespace MeshToolsLib
{
static void trackSurface(MeshLib::Element const* 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);
            }
        }
    }
}
 
bool MeshValidation::allNodesUsed(MeshLib::Mesh const& mesh)
{
    INFO("Looking for unused nodes...");
    MeshLib::NodeSearch ns(mesh);
    ns.searchUnused();
    if (!ns.getSearchedNodeIDs().empty())
    {
        INFO("{:d} unused mesh nodes found.", ns.getSearchedNodeIDs().size());
        return false;
    }
    return true;
}
 
bool MeshValidation::existCollapsibleNodes(MeshLib::Mesh& mesh,
                                           double const eps)
{
    MeshRevision const rev(mesh);
    INFO("Found {:d} potentially collapsible nodes.",
         rev.getNumberOfCollapsibleNodes(eps));
    return (rev.getNumberOfCollapsibleNodes() > 0);
}
 
void MeshValidation::evaluateElementGeometry(MeshLib::Mesh const& mesh)
{
    const std::vector<ElementErrorCode> codes(
        MeshToolsLib::MeshValidation::testElementGeometry(mesh));
    std::array<std::string,
               static_cast<std::size_t>(ElementErrorFlag::MaxValue)>
        output_str(MeshToolsLib::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 const 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));
    const 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 =
        MeshToolsLib::BoundaryExtraction::getBoundaryElementsAsMesh(
            mesh,
            MeshLib::getBulkIDString(MeshLib::MeshItemType::Node),
            MeshLib::getBulkIDString(MeshLib::MeshItemType::Cell),
            MeshLib::getBulkIDString(MeshLib::MeshItemType::Face));
    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);
}
}  // namespace MeshToolsLib