/** * \file * \author Karsten Rink * \date 2012-05-02 * \brief Implementation of the Element 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 "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) { if (e == this || e == nullptr || e->getDimension() != this->getDimension()) { 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); const unsigned dim(this->getDimension()); 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) { _neighbors[this->identifyFace(face_nodes)] = e; return std::optional<unsigned>(e->identifyFace(face_nodes)); } } } } return std::optional<unsigned>(); } bool Element::isBoundaryElement() const { return std::any_of( _neighbors, _neighbors + this->getNumberOfNeighbors(), [](MeshLib::Element const* const e) { return e == nullptr; }); } #ifndef NDEBUG 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: " << static_cast<std::underlying_type<MeshElemType>::type>(e.getGeomType()) << " with " << nnodes << " nodes: { "; for (unsigned n = 0; n < nnodes; ++n) { os << nodes[n]->getID() << " @ " << nodes[n] << " "; } os << "}\n"; return os; } #endif // NDEBUG 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(); } MeshLib::Node getCenterOfGravity(Element const& element) { const unsigned nNodes(element.getNumberOfBaseNodes()); MeshLib::Node center(0, 0, 0); for (unsigned i = 0; i < nNodes; ++i) { center[0] += (*element.getNode(i))[0]; center[1] += (*element.getNode(i))[1]; center[2] += (*element.getNode(i))[2]; } center[0] /= nNodes; center[1] /= nNodes; center[2] /= 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