LIE/Common/MeshUtils.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include "MeshUtils.h"
 
#include "BaseLib/Algorithm.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshSearch/NodeSearch.h"
#include "MeshLib/Node.h"
 
namespace ProcessLib
{
namespace LIE
{
namespace
{
// A class to check whether a node is located on a crack tip with
// the following conditions:
// - the number of connected fracture elements is one
// - the node is not located on a domain boundary
class IsCrackTip
{
public:
    explicit IsCrackTip(MeshLib::Mesh const& mesh)
        : _mesh(mesh), _fracture_element_dim(_mesh.getDimension() - 1)
    {
        _is_internal_node.resize(_mesh.getNumberOfNodes(), true);
 
        MeshLib::NodeSearch nodeSearch(mesh);
        nodeSearch.searchBoundaryNodes();
        for (auto i : nodeSearch.getSearchedNodeIDs())
        {
            _is_internal_node[i] = false;
        }
    }
 
    bool operator()(MeshLib::Node const& node) const
    {
        if (!_is_internal_node[node.getID()] ||
            !isBaseNode(node, _mesh.getElementsConnectedToNode(node)))
        {
            return false;
        }
 
        auto const elements_connected_to_node =
            _mesh.getElementsConnectedToNode(node);
        auto const n_connected_fracture_elements =
            count_if(cbegin(elements_connected_to_node),
                     cend(elements_connected_to_node),
                     [&](auto* const e)
                     { return e->getDimension() == _fracture_element_dim; });
        assert(n_connected_fracture_elements > 0);
 
        return (n_connected_fracture_elements == 1);
    }
 
private:
    MeshLib::Mesh const& _mesh;
    unsigned const _fracture_element_dim;
    std::vector<bool> _is_internal_node;
};
 
void findFracutreIntersections(
    MeshLib::Mesh const& mesh, std::vector<int> const& vec_fracture_mat_IDs,
    std::vector<std::vector<MeshLib::Node*>> const& vec_fracture_nodes,
    std::vector<std::vector<MeshLib::Element*>>& intersected_fracture_elements,
    std::vector<std::pair<std::size_t, std::vector<int>>>&
        vec_branch_nodeID_matIDs,
    std::vector<std::pair<std::size_t, std::vector<int>>>&
        vec_junction_nodeID_matIDs)
{
    auto const n_fractures = vec_fracture_mat_IDs.size();
    std::map<unsigned, unsigned> matID_to_fid;
    for (unsigned i = 0; i < n_fractures; i++)
    {
        matID_to_fid[vec_fracture_mat_IDs[i]] = i;
    }
 
    // make a vector all fracture nodes
    std::vector<std::size_t> all_fracture_nodes;
    for (auto& vec : vec_fracture_nodes)
    {
        transform(cbegin(vec), cend(vec), back_inserter(all_fracture_nodes),
                  [](auto* const n) { return n->getID(); });
    }
 
    // create a table of a node id and connected material IDs
    std::map<std::size_t, std::vector<std::size_t>> frac_nodeID_to_matIDs;
    for (unsigned i = 0; i < n_fractures; i++)
    {
        for (auto* node : vec_fracture_nodes[i])
        {
            frac_nodeID_to_matIDs[node->getID()].push_back(
                vec_fracture_mat_IDs[i]);
        }
    }
 
    auto const* const opt_material_ids = MeshLib::materialIDs(mesh);
 
    // find branch/junction nodes which connect to multiple fractures
    intersected_fracture_elements.resize(n_fractures);
    for (auto frac_nodeID_to_matID : frac_nodeID_to_matIDs)
    {
        auto nodeID = frac_nodeID_to_matID.first;
        auto const* node = mesh.getNode(frac_nodeID_to_matID.first);
        auto const& matIDs = frac_nodeID_to_matID.second;
        if (matIDs.size() < 2)
        {
            continue;  // no intersection
        }
 
        auto const elements_connected_to_node =
            mesh.getElementsConnectedToNode(*node);
        std::vector<MeshLib::Element*> conn_fracture_elements;
        {
            for (auto const* e : elements_connected_to_node)
            {
                if (e->getDimension() == (mesh.getDimension() - 1))
                {
                    conn_fracture_elements.push_back(
                        const_cast<MeshLib::Element*>(e));
                }
            }
        }
 
        std::map<int, int> vec_matID_counts;
        {
            for (auto matid : matIDs)
            {
                vec_matID_counts[matid] = 0;
            }
 
            for (auto const* e : conn_fracture_elements)
            {
                auto matid = (*opt_material_ids)[e->getID()];
                vec_matID_counts[matid]++;
            }
        }
 
        for (auto matid : matIDs)
        {
            auto fid = matID_to_fid[matid];
            for (auto* e : conn_fracture_elements)
            {
                auto e_matid = (*opt_material_ids)[e->getID()];
                if (matID_to_fid[e_matid] != fid)
                {  // only slave elements
                    intersected_fracture_elements[fid].push_back(e);
                }
            }
        }
 
        bool isBranch = false;
        {
            for (auto vec_matID_count : vec_matID_counts)
            {
                auto count = vec_matID_count.second;
                if (count % 2 == 1)
                {
                    isBranch = true;
                    break;
                }
            }
        }
 
        if (isBranch)
        {
            std::vector<int> branch_matIDs(2);
            for (auto vec_matID_count : vec_matID_counts)
            {
                auto matid = vec_matID_count.first;
                auto count = vec_matID_count.second;
                if (count % 2 == 0)
                {
                    branch_matIDs[0] = matid;  // master
                }
                else
                {
                    branch_matIDs[1] = matid;  // slave
                }
            }
            vec_branch_nodeID_matIDs.emplace_back(nodeID, branch_matIDs);
        }
        else
        {
            std::vector<int> junction_matIDs(2);
            junction_matIDs[0] = std::min(vec_matID_counts.begin()->first,
                                          vec_matID_counts.rbegin()->first);
            junction_matIDs[1] = std::max(vec_matID_counts.begin()->first,
                                          vec_matID_counts.rbegin()->first);
            vec_junction_nodeID_matIDs.emplace_back(nodeID, junction_matIDs);
        }
    }
 
    for (auto& elements : intersected_fracture_elements)
    {
        BaseLib::makeVectorUnique(elements);
    }
 
    DBUG("-> found {:d} branches and {:d} junctions",
         vec_branch_nodeID_matIDs.size(), vec_junction_nodeID_matIDs.size());
}
 
}  // namespace
 
void getFractureMatrixDataInMesh(
    MeshLib::Mesh const& mesh,
    std::vector<MeshLib::Element*>& vec_matrix_elements,
    std::vector<int>& vec_fracture_mat_IDs,
    std::vector<std::vector<MeshLib::Element*>>& vec_fracture_elements,
    std::vector<std::vector<MeshLib::Element*>>& vec_fracture_matrix_elements,
    std::vector<std::vector<MeshLib::Node*>>& vec_fracture_nodes,
    std::vector<std::pair<std::size_t, std::vector<int>>>&
        vec_branch_nodeID_matIDs,
    std::vector<std::pair<std::size_t, std::vector<int>>>&
        vec_junction_nodeID_matIDs)
{
    IsCrackTip isCrackTip(mesh);
 
    // get vectors of matrix elements and fracture elements
    vec_matrix_elements.reserve(mesh.getNumberOfElements());
    std::vector<MeshLib::Element*> all_fracture_elements;
    for (MeshLib::Element* e : mesh.getElements())
    {
        if (e->getDimension() == mesh.getDimension())
        {
            vec_matrix_elements.push_back(e);
        }
        else
        {
            all_fracture_elements.push_back(e);
        }
    }
    DBUG("-> found total {:d} matrix elements and {:d} fracture elements",
         vec_matrix_elements.size(), all_fracture_elements.size());
 
    // get fracture material IDs
    auto const material_ids = materialIDs(mesh);
    if (!material_ids)
    {
        OGS_FATAL("Could not access MaterialIDs property from mesh.");
    }
    transform(cbegin(all_fracture_elements), cend(all_fracture_elements),
              back_inserter(vec_fracture_mat_IDs),
              [&material_ids](auto const* e)
              { return (*material_ids)[e->getID()]; });
 
    BaseLib::makeVectorUnique(vec_fracture_mat_IDs);
    DBUG("-> found {:d} fracture material groups", vec_fracture_mat_IDs.size());
 
    // create a vector of fracture elements for each material
    vec_fracture_elements.resize(vec_fracture_mat_IDs.size());
    for (unsigned frac_id = 0; frac_id < vec_fracture_mat_IDs.size(); frac_id++)
    {
        const auto frac_mat_id = vec_fracture_mat_IDs[frac_id];
        std::vector<MeshLib::Element*>& vec_elements =
            vec_fracture_elements[frac_id];
        std::copy_if(all_fracture_elements.begin(), all_fracture_elements.end(),
                     std::back_inserter(vec_elements),
                     [&](MeshLib::Element const* const e)
                     { return (*material_ids)[e->getID()] == frac_mat_id; });
        DBUG("-> found {:d} elements on the fracture {:d}", vec_elements.size(),
             frac_id);
    }
 
    // get a vector of fracture nodes for each material
    vec_fracture_nodes.resize(vec_fracture_mat_IDs.size());
    for (unsigned frac_id = 0; frac_id < vec_fracture_mat_IDs.size(); frac_id++)
    {
        std::vector<MeshLib::Node*>& vec_nodes = vec_fracture_nodes[frac_id];
        for (MeshLib::Element const* const e : vec_fracture_elements[frac_id])
        {
            for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
            {
                if (isCrackTip(*e->getNode(i)))
                {
                    continue;
                }
                vec_nodes.push_back(const_cast<MeshLib::Node*>(e->getNode(i)));
            }
        }
        BaseLib::makeVectorUnique(vec_nodes,
                                  MeshLib::idsComparator<MeshLib::Node*>);
        DBUG("-> found {:d} nodes on the fracture {:d}", vec_nodes.size(),
             frac_id);
    }
 
    // find branch/junction nodes which connect to multiple fractures
    std::vector<std::vector<MeshLib::Element*>> intersected_fracture_elements;
    findFracutreIntersections(mesh, vec_fracture_mat_IDs, vec_fracture_nodes,
                              intersected_fracture_elements,
                              vec_branch_nodeID_matIDs,
                              vec_junction_nodeID_matIDs);
 
    // create a vector fracture elements and connected matrix elements,
    // which are passed to a DoF table
    for (unsigned fid = 0; fid < vec_fracture_elements.size(); fid++)
    {
        auto const& fracture_elements = vec_fracture_elements[fid];
        std::vector<MeshLib::Element*> vec_ele;
        // first, collect matrix elements
        for (MeshLib::Element const* const e : fracture_elements)
        {
            // it is sufficient to iterate over base nodes, because they are
            // already connected to all neighbours
            for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
            {
                MeshLib::Node const* node = e->getNode(i);
                if (isCrackTip(*node))
                {
                    continue;
                }
                auto const& elements_connected_to_node =
                    mesh.getElementsConnectedToNode(*node);
                for (unsigned j = 0; j < elements_connected_to_node.size(); j++)
                {
                    // only matrix elements
                    if (elements_connected_to_node[j]->getDimension() <
                        mesh.getDimension())
                    {
                        continue;
                    }
                    vec_ele.push_back(const_cast<MeshLib::Element*>(
                        elements_connected_to_node[j]));
                }
            }
        }
        BaseLib::makeVectorUnique(vec_ele,
                                  MeshLib::idsComparator<MeshLib::Element*>);
 
        // second, append fracture elements
        std::copy(fracture_elements.begin(), fracture_elements.end(),
                  std::back_inserter(vec_ele));
        // thirdly, append intersected fracture elements
        std::copy(intersected_fracture_elements[fid].begin(),
                  intersected_fracture_elements[fid].end(),
                  std::back_inserter(vec_ele));
 
        vec_fracture_matrix_elements.push_back(vec_ele);
    }
}
 
}  // namespace LIE
}  // namespace ProcessLib