Element.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include "Element.h"
 
#include "BaseLib/Logging.h"
#include "Line.h"
#include "MathLib/GeometricBasics.h"
#include "MathLib/MathTools.h"
#include "MeshLib/Node.h"
 
namespace MeshLib
{
Element::Element(std::size_t id) : _id(id), _neighbors(nullptr) {}
 
Element::~Element()
{
    delete[] this->_neighbors;
}
 
void Element::setNeighbor(Element* neighbor, unsigned const face_id)
{
    if (neighbor == this)
    {
        return;
    }
 
    this->_neighbors[face_id] = neighbor;
}
 
std::optional<unsigned> Element::addNeighbor(Element* e)
{
    const unsigned dim(this->getDimension());
    if (e == this || e == nullptr || e->getDimension() != dim)
    {
        return std::optional<unsigned>();
    }
 
    if (areNeighbors(this, e))
    {
        return std::optional<unsigned>();
    }
 
    Node const* face_nodes[3];
    const unsigned nNodes(this->getNumberOfBaseNodes());
    const unsigned eNodes(e->getNumberOfBaseNodes());
    const Node* const* e_nodes = e->getNodes();
    unsigned count(0);
    for (unsigned i(0); i < nNodes; i++)
    {
        for (unsigned j(0); j < eNodes; j++)
        {
            if (getNode(i) == e_nodes[j])
            {
                face_nodes[count] = getNode(i);
                // increment shared nodes counter and check if enough nodes are
                // similar to be sure e is a neighbour of this
                if ((++count) >= dim)
                {
                    auto const this_face_id = this->identifyFace(face_nodes);
                    if (this_face_id == std::numeric_limits<unsigned>::max())
                    {
                        OGS_FATAL(
                            "Could not find face of {} for the face nodes [{}, "
                            "{}, {}].",
                            *this, *face_nodes[0], *face_nodes[1],
                            *face_nodes[2]);
                    }
                    _neighbors[this_face_id] = e;
                    auto const neighbors_face_id = e->identifyFace(face_nodes);
                    if (neighbors_face_id ==
                        std::numeric_limits<unsigned>::max())
                    {
                        OGS_FATAL(
                            "Could not find face of the neighbour {} for the "
                            "face nodes [{}, {}, {}].\n"
                            "For the {} the face id {} was identified.",
                            *e, *face_nodes[0], *face_nodes[1], *face_nodes[2],
                            *this, this_face_id);
                    }
                    return neighbors_face_id;
                }
            }
        }
    }
 
    return std::optional<unsigned>();
}
 
bool Element::isBoundaryElement() const
{
    return std::any_of(_neighbors, _neighbors + this->getNumberOfNeighbors(),
                       [](MeshLib::Element const* const e)
                       { return e == nullptr; });
}
 
std::ostream& operator<<(std::ostream& os, Element const& e)
{
    os << "element #" << e._id << " @ " << &e << " with "
       << e.getNumberOfNeighbors() << " neighbours\n";
 
    unsigned const nnodes = e.getNumberOfNodes();
    MeshLib::Node* const* const nodes = e.getNodes();
    os << "MeshElemType: " << MeshLib::MeshElemType2String(e.getGeomType())
       << " with " << nnodes << " nodes: {\n";
    for (unsigned n = 0; n < nnodes; ++n)
    {
        os << "  " << *nodes[n] << " @ " << nodes[n] << '\n';
    }
    return os << '}';
}
 
bool areNeighbors(Element const* const element, Element const* const other)
{
    unsigned nNeighbors(element->getNumberOfNeighbors());
    for (unsigned i = 0; i < nNeighbors; i++)
    {
        if (element->getNeighbor(i) == other)
        {
            return true;
        }
    }
    return false;
}
 
bool hasZeroVolume(MeshLib::Element const& element)
{
    return element.getContent() < std::numeric_limits<double>::epsilon();
}
 
MathLib::Point3d getCenterOfGravity(Element const& element)
{
    const unsigned nNodes(element.getNumberOfBaseNodes());
    MathLib::Point3d center{{0, 0, 0}};
    for (unsigned i = 0; i < nNodes; ++i)
    {
        center.asEigenVector3d() += element.getNode(i)->asEigenVector3d();
    }
    center.asEigenVector3d() /= nNodes;
    return center;
}
 
std::pair<double, double> computeSqrNodeDistanceRange(
    MeshLib::Element const& element, bool const check_allnodes)
{
    double min = std::numeric_limits<double>::max();
    double max = 0;
    const unsigned nnodes = check_allnodes ? element.getNumberOfNodes()
                                           : element.getNumberOfBaseNodes();
    for (unsigned i = 0; i < nnodes; i++)
    {
        for (unsigned j = i + 1; j < nnodes; j++)
        {
            const double dist(
                MathLib::sqrDist(*element.getNode(i), *element.getNode(j)));
            min = std::min(dist, min);
            max = std::max(dist, max);
        }
    }
    return {min, max};
}
 
std::pair<double, double> computeSqrEdgeLengthRange(Element const& element)
{
    double min = std::numeric_limits<double>::max();
    double max = 0;
    const unsigned nEdges(element.getNumberOfEdges());
    for (unsigned i = 0; i < nEdges; i++)
    {
        const double dist(MathLib::sqrDist(*element.getEdgeNode(i, 0),
                                           *element.getEdgeNode(i, 1)));
        min = std::min(dist, min);
        max = std::max(dist, max);
    }
    return {min, max};
}
 
bool isPointInElementXY(MathLib::Point3d const& p, Element const& e)
{
    for (std::size_t i(0); i < e.getNumberOfBaseNodes(); ++i)
    {
        if (MathLib::sqrDist2d(p, *e.getNode(i)) <
            std::numeric_limits<double>::epsilon())
        {
            return true;
        }
    }
 
    if (e.getGeomType() == MeshElemType::TRIANGLE)
    {
        MathLib::Point3d const& n0(*e.getNode(0));
        MathLib::Point3d const& n1(*e.getNode(1));
        MathLib::Point3d const& n2(*e.getNode(2));
 
        return MathLib::isPointInTriangleXY(p, n0, n1, n2);
    }
    if (e.getGeomType() == MeshElemType::QUAD)
    {
        MathLib::Point3d const& n0(*e.getNode(0));
        MathLib::Point3d const& n1(*e.getNode(1));
        MathLib::Point3d const& n2(*e.getNode(2));
        MathLib::Point3d const& n3(*e.getNode(3));
 
        return MathLib::isPointInTriangleXY(p, n0, n1, n2) ||
               MathLib::isPointInTriangleXY(p, n0, n2, n3);
    }
 
    WARN("isPointInElementXY: element type '{:s}' is not supported.",
         MeshLib::MeshElemType2String(e.getGeomType()));
    return false;
}
 
unsigned getNodeIDinElement(Element const& element, const MeshLib::Node* node)
{
    const unsigned nNodes(element.getNumberOfNodes());
    for (unsigned i(0); i < nNodes; i++)
    {
        if (node == element.getNode(i))
        {
            return i;
        }
    }
    return std::numeric_limits<unsigned>::max();
}
 
std::size_t getNodeIndex(Element const& element, unsigned const idx)
{
#ifndef NDEBUG
    if (idx >= element.getNumberOfNodes())
    {
        ERR("Error in MeshLib::getNodeIndex() - Index does not "
            "exist.");
        return std::numeric_limits<std::size_t>::max();
    }
#endif
    return element.getNode(idx)->getID();
}
 
}  // namespace MeshLib