/** * \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 "MathLib/GeometricBasics.h" #include "MathLib/MathTools.h" #include "MeshLib/Node.h" #include "Line.h" namespace MeshLib { Element::Element(std::size_t id) : _nodes(nullptr), _id(id), _content(-1.0), _neighbors(nullptr) { } Element::~Element() { delete [] this->_nodes; delete [] this->_neighbors; } void Element::setNeighbor(Element* neighbor, unsigned const face_id) { if (neighbor == this) { return; } this->_neighbors[face_id] = neighbor; } boost::optional Element::addNeighbor(Element* e) { if (e == this || e == nullptr || e->getDimension() != this->getDimension()) { return boost::optional(); } if (areNeighbors(this, e)) { return boost::optional(); } Node* 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 (_nodes[i] == e_nodes[j]) { face_nodes[count] = _nodes[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 boost::optional(e->identifyFace(face_nodes)); } } } } return boost::optional(); } const Element* Element::getNeighbor(unsigned i) const { #ifndef NDEBUG if (i < getNumberOfNeighbors()) #endif { return _neighbors[i]; } #ifndef NDEBUG ERR("Error in MeshLib::Element::getNeighbor() - Index does not exist."); return nullptr; #endif } unsigned Element::getNodeIDinElement(const MeshLib::Node* node) const { const unsigned nNodes (this->getNumberOfNodes()); for (unsigned i(0); i < nNodes; i++) { if (node == _nodes[i]) { return i; } } return std::numeric_limits::max(); } const Node* Element::getNode(unsigned i) const { #ifndef NDEBUG if (i < getNumberOfNodes()) #endif { return _nodes[i]; } #ifndef NDEBUG ERR("Error in MeshLib::Element::getNode() - Index {:d} in {:s}", i, MeshElemType2String(getGeomType())); return nullptr; #endif } void Element::setNode(unsigned idx, Node* node) { #ifndef NDEBUG if (idx < getNumberOfNodes()) #endif { _nodes[idx] = node; } } std::size_t Element::getNodeIndex(unsigned i) const { #ifndef NDEBUG if (i < getNumberOfNodes()) #endif { return _nodes[i]->getID(); } #ifndef NDEBUG ERR("Error in MeshLib::Element::getNodeIndex() - Index does not exist."); return std::numeric_limits::max(); #endif } 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::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::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 computeSqrNodeDistanceRange( MeshLib::Element const& element, bool const check_allnodes) { double min = std::numeric_limits::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 computeSqrEdgeLengthRange(Element const& element) { double min = std::numeric_limits::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::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; } } // namespace MeshLib