Merge branch 'RefactorMeshLibElementRules' into 'master'

Refactor MeshLib element rules

See merge request ogs/ogs!4675

MeshLib/Elements/CellRule.cpp

@@ -12,7 +12,6 @@
 
 #include "Element.h"
 #include "FaceRule.h"
-#include "MeshLib/Node.h"
 
 namespace MeshLib
 {

MeshLib/Elements/CellRule.h

@@ -10,6 +10,8 @@
 
 #pragma once
 
+#include "MeshLib/Node.h"
+
 namespace MeshLib
 {
 
@@ -30,6 +32,37 @@ class CellRule
  * center of gravity to lie outside of the actual element
  */
  static bool testElementNodeOrder(Element const& e);
+
+protected:
+ /// Returns the ID of a face given an array of nodes.
+ template <typename ElementRule>
+ static unsigned identifyFace(Node const* const* element_nodes,
+ Node const* nodes[ElementRule::dimension])
+ {
+ for (unsigned i = 0; i < ElementRule::n_faces; i++)
+ {
+ unsigned flag(0);
+ constexpr std::size_t n = sizeof(ElementRule::face_nodes[0]) /
+ sizeof(ElementRule::face_nodes[0][0]);
+ for (unsigned j = 0; j < n; j++)
+ {
+ for (unsigned k = 0; k < ElementRule::dimension; k++)
+ {
+ if (ElementRule::face_nodes[i][j] != 99 &&
+ element_nodes[ElementRule::face_nodes[i][j]] ==
+ nodes[k])
+ {
+ flag++;
+ }
+ }
+ }
+ if (flag == ElementRule::dimension)
+ {
+ return i;
+ }
+ }
+ return std::numeric_limits<unsigned>::max();
+ }
 }; /* class */
 
 } // namespace MeshLib

MeshLib/Elements/FaceRule.h

@@ -42,6 +42,36 @@ class FaceRule
  /// Returns the surface normal of a 2D element.
  static Eigen::Vector3d getSurfaceNormal(Element const& e);
 
+protected:
+ /// Returns the ID of an edge given an array of nodes.
+ template <typename ElementRule>
+ static unsigned identifyFace(Node const* const* element_nodes,
+ Node const* nodes[ElementRule::dimension])
+ {
+ for (unsigned i = 0; i < ElementRule::n_edges; i++)
+ {
+ unsigned flag(0);
+ constexpr std::size_t n = sizeof(ElementRule::edge_nodes[0]) /
+ sizeof(ElementRule::edge_nodes[0][0]);
+ for (unsigned j = 0; j < n; j++)
+ {
+ for (unsigned k = 0; k < ElementRule::dimension; k++)
+ {
+ if (element_nodes[ElementRule::edge_nodes[i][j]] ==
+ nodes[k])
+ {
+ flag++;
+ }
+ }
+ }
+ if (flag == ElementRule::dimension)
+ {
+ return i;
+ }
+ }
+ return std::numeric_limits<unsigned>::max();
+ }
+
 }; /* class */
 
 } // namespace MeshLib

MeshLib/Elements/HexRule.cpp

@@ -0,0 +1,84 @@
+/**
+ * \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* element_nodes)
+{
+ return MathLib::calcTetrahedronVolume(*element_nodes[4],
+ *element_nodes[7],
+ *element_nodes[5],
+ *element_nodes[0]) +
+ MathLib::calcTetrahedronVolume(*element_nodes[5],
+ *element_nodes[3],
+ *element_nodes[1],
+ *element_nodes[0]) +
+ MathLib::calcTetrahedronVolume(*element_nodes[5],
+ *element_nodes[7],
+ *element_nodes[3],
+ *element_nodes[0]) +
+ MathLib::calcTetrahedronVolume(*element_nodes[5],
+ *element_nodes[7],
+ *element_nodes[6],
+ *element_nodes[2]) +
+ MathLib::calcTetrahedronVolume(*element_nodes[1],
+ *element_nodes[3],
+ *element_nodes[5],
+ *element_nodes[2]) +
+ MathLib::calcTetrahedronVolume(*element_nodes[3],
+ *element_nodes[7],
+ *element_nodes[5],
+ *element_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

MeshLib/Elements/HexRule.h

@@ -0,0 +1,53 @@
+/**
+ * \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* element_nodes);
+};
+} // namespace MeshLib

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 @@ class HexRule20 : public HexRule8
  /// 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* element_nodes,
+ Node const* nodes[3])
+ {
+ return CellRule::identifyFace<HexRule20>(element_nodes, nodes);
+ }
 }; /* class */
 
 } // namespace MeshLib

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