diff --git a/MeshLib/Elements/HexRule.cpp b/MeshLib/Elements/HexRule.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6210426f4a9a9b55ebbdfe2652d4af1c37c4ca9f --- /dev/null +++ b/MeshLib/Elements/HexRule.cpp @@ -0,0 +1,72 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, 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 "HexRule.h" + +#include "MathLib/GeometricBasics.h" +#include "MeshLib/Node.h" +#include "Quad.h" + +namespace MeshLib +{ +double HexRule::computeVolume(Node const* const* _nodes) +{ + return MathLib::calcTetrahedronVolume( + *_nodes[4], *_nodes[7], *_nodes[5], *_nodes[0]) + + MathLib::calcTetrahedronVolume( + *_nodes[5], *_nodes[3], *_nodes[1], *_nodes[0]) + + MathLib::calcTetrahedronVolume( + *_nodes[5], *_nodes[7], *_nodes[3], *_nodes[0]) + + MathLib::calcTetrahedronVolume( + *_nodes[5], *_nodes[7], *_nodes[6], *_nodes[2]) + + MathLib::calcTetrahedronVolume( + *_nodes[1], *_nodes[3], *_nodes[5], *_nodes[2]) + + MathLib::calcTetrahedronVolume( + *_nodes[3], *_nodes[7], *_nodes[5], *_nodes[2]); +} + +bool HexRule::isPntInElement(Node const* const* nodes, + MathLib::Point3d const& pnt, + double eps) +{ + return (MathLib::isPointInTetrahedron( + pnt, *nodes[4], *nodes[7], *nodes[5], *nodes[0], eps) || + MathLib::isPointInTetrahedron( + pnt, *nodes[5], *nodes[3], *nodes[1], *nodes[0], eps) || + MathLib::isPointInTetrahedron( + pnt, *nodes[5], *nodes[7], *nodes[3], *nodes[0], eps) || + MathLib::isPointInTetrahedron( + pnt, *nodes[5], *nodes[7], *nodes[6], *nodes[2], eps) || + MathLib::isPointInTetrahedron( + pnt, *nodes[1], *nodes[3], *nodes[5], *nodes[2], eps) || + MathLib::isPointInTetrahedron( + pnt, *nodes[3], *nodes[7], *nodes[5], *nodes[2], eps)); +} + +ElementErrorCode HexRule::validate(const Element* e) +{ + ElementErrorCode error_code; + error_code[ElementErrorFlag::ZeroVolume] = hasZeroVolume(*e); + + for (unsigned i = 0; i < 6; ++i) + { + if (error_code.all()) + { + break; + } + + const MeshLib::Element* quad(e->getFace(i)); + error_code |= quad->validate(); + delete quad; + } + error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder(); + return error_code; +} + +} // end namespace MeshLib diff --git a/MeshLib/Elements/HexRule.h b/MeshLib/Elements/HexRule.h new file mode 100644 index 0000000000000000000000000000000000000000..203c3b9c3a62c4a04e094ea858c6aed348a76ceb --- /dev/null +++ b/MeshLib/Elements/HexRule.h @@ -0,0 +1,80 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#pragma once + +#include "CellRule.h" +#include "Element.h" +#include "MeshLib/MeshEnums.h" + +namespace MeshLib +{ +class HexRule : public CellRule +{ +public: + /// Constant: The number of base nodes for this element + static const unsigned n_base_nodes = 8u; + + /// Constant: The geometric type of the element + static const MeshElemType mesh_elem_type = MeshElemType::HEXAHEDRON; + + /// Constant: The number of faces + static const unsigned n_faces = 6; + + /// Constant: The number of edges + static const unsigned n_edges = 12; + + /// Constant: The number of neighbors + static const unsigned n_neighbors = 6; + + /** + * \copydoc MeshLib::Element::isPntInElement() + * \param nodes the nodes of the element. + */ + static bool isPntInElement(Node const* const* nodes, + MathLib::Point3d const& pnt, double eps); + + /** + * Tests if the element is geometrically valid. + */ + static ElementErrorCode validate(const Element* e); + + /// Calculates the volume of a convex hexahedron by partitioning it into six + /// tetrahedra. + static double computeVolume(Node const* const* _nodes); + +protected: + template <std::size_t N> + static unsigned identifyFace(Node const* const* _nodes, + Node const* nodes[3], + unsigned const face_nodes[6][N]) + { + for (unsigned i = 0; i < 6; i++) + { + unsigned flag(0); + for (unsigned j = 0; j < 4; j++) + { + for (unsigned k = 0; k < 3; k++) + { + if (_nodes[face_nodes[i][j]] == nodes[k]) + { + flag++; + } + } + } + if (flag == 3) + { + return i; + } + } + return std::numeric_limits<unsigned>::max(); + } +}; +} // namespace MeshLib diff --git a/MeshLib/Elements/HexRule20.h b/MeshLib/Elements/HexRule20.h index 19bb18bacb9a68a9950825d30a2a88f3e3d21362..b5024f5aaa976050e7cd88d903b26e64ea6fb014 100644 --- a/MeshLib/Elements/HexRule20.h +++ b/MeshLib/Elements/HexRule20.h @@ -10,10 +10,8 @@ #pragma once -#include "MeshLib/MeshEnums.h" -#include "Element.h" #include "EdgeReturn.h" -#include "HexRule8.h" +#include "HexRule.h" namespace MeshLib { @@ -43,7 +41,7 @@ namespace MeshLib * * \endcode */ -class HexRule20 : public HexRule8 +class HexRule20 : public HexRule { public: /// Constant: The number of all nodes for this element @@ -64,6 +62,12 @@ public: /// Returns the i-th face of the element. static const Element* getFace(const Element* e, unsigned i); + /// Returns the ID of a face given an array of nodes. + static unsigned identifyFace(Node const* const* _nodes, + Node const* nodes[3]) + { + return HexRule::identifyFace(_nodes, nodes, face_nodes); + } }; /* class */ } // namespace MeshLib diff --git a/MeshLib/Elements/HexRule8.cpp b/MeshLib/Elements/HexRule8.cpp index 1144c14a71e04639552e5b67f5ca7e0f74aa2405..52d31ad4d1e31a8c54d253ac58b19d7a673e6826 100644 --- a/MeshLib/Elements/HexRule8.cpp +++ b/MeshLib/Elements/HexRule8.cpp @@ -58,82 +58,4 @@ const Element* HexRule8::getFace(const Element* e, unsigned i) ERR("Error in MeshLib::Element::getFace() - Index {:d} does not exist.", i); return nullptr; } - -double HexRule8::computeVolume(Node const* const* _nodes) -{ - return MathLib::calcTetrahedronVolume( - *_nodes[4], *_nodes[7], *_nodes[5], *_nodes[0]) + - MathLib::calcTetrahedronVolume( - *_nodes[5], *_nodes[3], *_nodes[1], *_nodes[0]) + - MathLib::calcTetrahedronVolume( - *_nodes[5], *_nodes[7], *_nodes[3], *_nodes[0]) + - MathLib::calcTetrahedronVolume( - *_nodes[5], *_nodes[7], *_nodes[6], *_nodes[2]) + - MathLib::calcTetrahedronVolume( - *_nodes[1], *_nodes[3], *_nodes[5], *_nodes[2]) + - MathLib::calcTetrahedronVolume( - *_nodes[3], *_nodes[7], *_nodes[5], *_nodes[2]); -} - -bool HexRule8::isPntInElement(Node const* const* nodes, - MathLib::Point3d const& pnt, - double eps) -{ - return (MathLib::isPointInTetrahedron( - pnt, *nodes[4], *nodes[7], *nodes[5], *nodes[0], eps) || - MathLib::isPointInTetrahedron( - pnt, *nodes[5], *nodes[3], *nodes[1], *nodes[0], eps) || - MathLib::isPointInTetrahedron( - pnt, *nodes[5], *nodes[7], *nodes[3], *nodes[0], eps) || - MathLib::isPointInTetrahedron( - pnt, *nodes[5], *nodes[7], *nodes[6], *nodes[2], eps) || - MathLib::isPointInTetrahedron( - pnt, *nodes[1], *nodes[3], *nodes[5], *nodes[2], eps) || - MathLib::isPointInTetrahedron( - pnt, *nodes[3], *nodes[7], *nodes[5], *nodes[2], eps)); -} - -unsigned HexRule8::identifyFace(Node const* const* _nodes, Node const* nodes[3]) -{ - for (unsigned i = 0; i < 6; i++) - { - unsigned flag(0); - for (unsigned j = 0; j < 4; j++) - { - for (unsigned k = 0; k < 3; k++) - { - if (_nodes[face_nodes[i][j]] == nodes[k]) - { - flag++; - } - } - } - if (flag == 3) - { - return i; - } - } - return std::numeric_limits<unsigned>::max(); -} - -ElementErrorCode HexRule8::validate(const Element* e) -{ - ElementErrorCode error_code; - error_code[ElementErrorFlag::ZeroVolume] = hasZeroVolume(*e); - - for (unsigned i = 0; i < 6; ++i) - { - if (error_code.all()) - { - break; - } - - const MeshLib::Element* quad(e->getFace(i)); - error_code |= quad->validate(); - delete quad; - } - error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder(); - return error_code; -} - } // end namespace MeshLib diff --git a/MeshLib/Elements/HexRule8.h b/MeshLib/Elements/HexRule8.h index d1501fbfb493bf721591370fe0f4d1f579f783c7..d1bb1ef2c79d9f43e1a4a4a540678c94f01a8d10 100644 --- a/MeshLib/Elements/HexRule8.h +++ b/MeshLib/Elements/HexRule8.h @@ -10,14 +10,11 @@ #pragma once -#include "MeshLib/MeshEnums.h" -#include "Element.h" #include "EdgeReturn.h" -#include "CellRule.h" +#include "HexRule.h" namespace MeshLib { - /** * A 8-nodes Hexahedron Element. * \code @@ -43,30 +40,15 @@ namespace MeshLib * * \endcode */ -class HexRule8 : public CellRule +class HexRule8 : public HexRule { public: - /// Constant: The number of base nodes for this element - static const unsigned n_base_nodes = 8u; - /// Constant: The number of all nodes for this element static const unsigned n_all_nodes = 8u; - /// Constant: The geometric type of the element - static const MeshElemType mesh_elem_type = MeshElemType::HEXAHEDRON; - /// Constant: The FEM type of the element static const CellType cell_type = CellType::HEX8; - /// Constant: The number of faces - static const unsigned n_faces = 6; - - /// Constant: The number of edges - static const unsigned n_edges = 12; - - /// Constant: The number of neighbors - static const unsigned n_neighbors = 6; - /// Constant: Local node index table for faces static const unsigned face_nodes[6][4]; @@ -79,25 +61,12 @@ public: /// Returns the i-th face of the element. static const Element* getFace(const Element* e, unsigned i); - /** - * \copydoc MeshLib::Element::isPntInElement() - * \param nodes the nodes of the element. - */ - static bool isPntInElement(Node const* const* nodes, - MathLib::Point3d const& pnt, double eps); - - /** - * Tests if the element is geometrically valid. - */ - static ElementErrorCode validate(const Element* e); - /// Returns the ID of a face given an array of nodes. - static unsigned identifyFace(Node const* const* /*_nodes*/, - Node const* nodes[3]); - - /// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra. - static double computeVolume(Node const* const* _nodes); - + static unsigned identifyFace(Node const* const* _nodes, + Node const* nodes[3]) + { + return HexRule::identifyFace(_nodes, nodes, face_nodes); + } }; /* class */ } // namespace MeshLib