feat: unify implementation of polyhedra

python/overlap/overlap.cpp

@@ -18,18 +18,17 @@
  * along with this program. If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include "overlap/overlap.hpp"
+
+#include <pybind11/eigen.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
 #include <cctype>
 #include <locale>
 #include <map>
 #include <typeindex>
-#include <typeinfo>
-
-#include "pybind11/eigen.h"
-#include "pybind11/numpy.h"
-#include "pybind11/pybind11.h"
-#include "pybind11/stl.h"
-
-#include "overlap/overlap.hpp"
 
 namespace py = pybind11;
 
@@ -42,64 +41,62 @@ using namespace overlap;
 
 template<typename Element,
          typename = std::enable_if_t<detail::is_element_v<Element>>>
-void createBindings(py::module& m) {
+void create_bindings(py::module& m) {
   static const auto element_names = std::map<std::type_index, std::string>{
       {std::type_index(typeid(Tetrahedron)), "Tetrahedron"},
       {std::type_index(typeid(Wedge)), "Wedge"},
       {std::type_index(typeid(Hexahedron)), "Hexahedron"}};
 
   const auto& name = element_names.at(std::type_index(typeid(Element)));
-  const auto nameLower =
+  const auto name_lower =
       std::tolower(name.front(), std::locale{}) + name.substr(1);
 
-  static constexpr auto nrVertices = detail::num_vertices<Element>();
+  static constexpr auto num_vertices = Element::num_vertices();
 
   py::class_<Element>(m, name.c_str())
-      .def(py::init<std::array<Vector, nrVertices>>(), py::arg("vertices"))
+      .def(py::init<std::array<Vector, num_vertices>>(), py::arg("vertices"))
       .def(py::init([](py::array_t<double> vertices) {
              auto proxy = vertices.unchecked<2>();
-             if (proxy.shape(0) != nrVertices || proxy.shape(1) != 3) {
+             if (proxy.shape(0) != num_vertices || proxy.shape(1) != 3) {
                throw std::invalid_argument{
                    "invalid shape for vertex list, must be (" +
-                   std::to_string(nrVertices) + ", 3)"};
+                   std::to_string(num_vertices) + ", 3)"};
              }
 
-             std::array<Vector, nrVertices> tmp{};
+             std::array<Vector, num_vertices> tmp{};
              for (py::ssize_t v = 0; v < proxy.shape(0); ++v) {
                tmp[v] = Vector{proxy(v, 0), proxy(v, 1), proxy(v, 2)};
              }
 
              return Element{std::move(tmp)};
            }),
            py::arg("vertices"))
-      .def_readonly("vertices", &Element::vertices,
-                    "Return the vertices of the element.")
-      .def_readonly("center", &Element::center,
-                    "Return the center point of the element.")
-      .def_readonly("volume", &Element::volume,
-                    "Return the volume of the element.")
-      .def_property_readonly(
-          "surface_area",
-          [](const Element& elem) { return elem.surface_area(); },
-          "Return the surface area of the element.");
+      .def_property_readonly("vertices", &Element::vertices,
+                             "Return the vertices of the element.")
+      .def_property_readonly("center", &Element::center,
+                             "Return the center point of the element.")
+      .def_property_readonly("volume", &Element::volume,
+                             "Return the volume of the element.")
+      .def_property_readonly("surface_area", &Element::surface_area,
+                             "Return the surface area of the element.");
 
   m.def(
       "overlap_volume",
       overload_cast_<const Sphere&, const Element&>()(&overlap_volume<Element>),
-      "sphere"_a, py::arg{nameLower.c_str()},
-      ("Calculate the overlap volume of a sphere and a " + nameLower + ".")
+      "sphere"_a, py::arg{name_lower.c_str()},
+      ("Calculate the overlap volume of a sphere and a " + name_lower + ".")
           .c_str());
 
   m.def("overlap_area",
         overload_cast_<const Sphere&, const Element&>()(&overlap_area<Element>),
-        "sphere"_a, py::arg{nameLower.c_str()},
-        ("Calculate the overlap area of a sphere and a " + nameLower + ".")
+        "sphere"_a, py::arg{name_lower.c_str()},
+        ("Calculate the overlap area of a sphere and a " + name_lower + ".")
             .c_str());
 }
 
 PYBIND11_MODULE(overlap, m) {
   m.doc() = R"pbdoc(
-        Pybind11 example plugin
+        Overlap Python bindings
         -----------------------
 
         This originates fom CPP
@@ -128,7 +125,7 @@ PYBIND11_MODULE(overlap, m) {
           "surface_area", [](const Sphere& s) { return s.surface_area(); },
           "Return the surface area of the sphere.");
 
-  createBindings<Tetrahedron>(m);
-  createBindings<Wedge>(m);
-  createBindings<Hexahedron>(m);
+  create_bindings<Tetrahedron>(m);
+  create_bindings<Wedge>(m);
+  create_bindings<Hexahedron>(m);
 }

tests/python/test_overlap.py

@@ -113,7 +113,7 @@ def test_wedge(self):
 
         np.testing.assert_almost_equal(overlap.overlap_area(sphere, wedge), [
             sphere.surface_area / 2,
-            0., 0., 0., sphere.radius**2 * np.pi, 0.,
+            0., 0., sphere.radius**2 * np.pi, 0., 0.,
             sphere.radius**2 * np.pi
         ])
 

tests/src/CMakeLists.txt

@@ -74,9 +74,9 @@ set(unit_tests
     clamp
     contains
     decompose_elements
-    elements
     normal_newell
     polygon
+    polyhedron
     regularized_wedge
     regularized_wedge_area
     sphere
@@ -107,11 +107,11 @@ if(_have_fenv_h)
 endif()
 cmake_reset_check_state()
 
-# activate additional compiler warnings
+# activate additional compiler warnings and keep frame pointers
 if(CMAKE_CXX_COMPILER_ID MATCHES "(GNU|Clang)" AND CMAKE_BUILD_TYPE MATCHES
                                                    "(Debug|RelWithDebInfo)"
 )
-  string(APPEND CMAKE_CXX_FLAGS " -Wall -Werror -Wextra")
+  string(APPEND CMAKE_CXX_FLAGS "-Wall -Werror -Wextra -fno-omit-frame-pointer")
 endif()
 
 # register the individual unit tests

tests/src/common.hpp

@@ -21,21 +21,19 @@
 #ifndef OVERLAP_TEST_COMMON_HPP
 #define OVERLAP_TEST_COMMON_HPP
 
-#include <iostream>
 #include <string_view>
 
 #include <doctest/doctest.h>
 
 class AssertionError : public std::runtime_error {
  public:
-  static inline constexpr auto assertion_failed_msg =
-      "[overlap] assertion failed: ";
+  static constexpr auto assertion_failed_msg = "[overlap] assertion failed: ";
 
   explicit AssertionError(std::string_view msg) :
       std::runtime_error{assertion_failed_msg + std::string{msg}} {}
 };
 
-[[noreturn]] inline void throw_assertion_error(std::string_view msg) {
+[[noreturn]] void throw_assertion_error(std::string_view msg) {
   throw AssertionError{msg};
 }
 
@@ -102,12 +100,6 @@ inline void validate_overlap_volume(const Sphere& s, const Hexahedron& hex,
     overlapCalcWedges += overlap_volume(s, wedge);
   }
 
-  std::cout << "volume hex:    " << overlapCalcHex << std::endl;
-  std::cout << "volume wedges: " << overlapCalcWedges << std::endl;
-  std::cout << "volume tets5:  " << overlapCalcTets5 << std::endl;
-  std::cout << "volume tets6:  " << overlapCalcTets6 << std::endl;
-  std::cout << "volume tets24: " << overlapCalcTets24 << std::endl;
-
   CHECK(overlapCalcHex == Approx(overlapCalcWedges).epsilon(epsilon));
   CHECK(overlapCalcHex == Approx(overlapCalcTets5).epsilon(epsilon));
   CHECK(overlapCalcHex == Approx(overlapCalcTets6).epsilon(epsilon));
@@ -129,15 +121,13 @@ inline void validate_overlap_area(const Sphere& s, const Hexahedron& hex,
 
   const auto areaCalcHex = overlap_area(s, hex)[0];
 
-  auto areaCalcWedges = Scalar{0};
   auto areaCalcTets5 = Scalar{0};
-  auto areaCalcTets6 = Scalar{0};
-  auto areaCalcTets24 = Scalar{0};
-
   for (const auto& tet : tets5) {
     areaCalcTets5 += overlap_area(s, tet)[0];
   }
 
+  auto areaCalcTets6 = Scalar{0};
+  auto areaCalcTets24 = Scalar{0};
   for (const auto& tet : tets6) {
     std::array<Tetrahedron, 4> subTets;
     decompose(tet, subTets);
@@ -149,19 +139,11 @@ inline void validate_overlap_area(const Sphere& s, const Hexahedron& hex,
     areaCalcTets6 += overlap_area(s, tet)[0];
   }
 
+  auto areaCalcWedges = Scalar{0};
   for (const auto& wedge : wedges) {
     areaCalcWedges += overlap_area(s, wedge)[0];
   }
 
-  std::cout << "sphere center: [" << s.center.transpose()
-            << "], radius: " << s.radius << std::endl;
-
-  std::cout << "sphere surface hex:    " << areaCalcHex << std::endl;
-  std::cout << "sphere surface wedges: " << areaCalcWedges << std::endl;
-  std::cout << "sphere surface tets5:  " << areaCalcTets5 << std::endl;
-  std::cout << "sphere surface tets6:  " << areaCalcTets6 << std::endl;
-  std::cout << "sphere surface tets24: " << areaCalcTets24 << std::endl;
-
   CHECK(areaCalcHex == Approx(areaCalcWedges).epsilon(epsilon));
   CHECK(areaCalcHex == Approx(areaCalcTets5).epsilon(epsilon));
   CHECK(areaCalcHex == Approx(areaCalcTets6).epsilon(epsilon));

tests/src/test_decompose_elements.cpp

@@ -35,28 +35,28 @@ TEST_SUITE("DecomposeElements") {
 
     auto tets5Volume = Scalar{0};
     for (const auto& tet : tets5) {
-      tets5Volume += tet.volume;
+      tets5Volume += tet.volume();
     }
 
     auto tets6Volume = Scalar{0};
     auto tets24Volume = Scalar{0};
     for (const auto& tet : tets6) {
       std::array<Tetrahedron, 4> subTets;
       decompose(tet, subTets);
-      tets6Volume += tet.volume;
+      tets6Volume += tet.volume();
 
       for (const auto& subTet : subTets) {
-        tets24Volume += subTet.volume;
+        tets24Volume += subTet.volume();
       }
     }
 
     constexpr auto epsilon =
         Scalar{5e2} * std::numeric_limits<Scalar>::epsilon();
 
-    REQUIRE(hex.volume == Approx(tets5Volume).epsilon(epsilon));
-    REQUIRE(hex.volume == Approx(tets6Volume).epsilon(epsilon));
-    REQUIRE(hex.volume == Approx(tets24Volume).epsilon(epsilon));
-    REQUIRE(hex.volume ==
-            Approx(wedges[0].volume + wedges[1].volume).epsilon(epsilon));
+    REQUIRE(hex.volume() == Approx(tets5Volume).epsilon(epsilon));
+    REQUIRE(hex.volume() == Approx(tets6Volume).epsilon(epsilon));
+    REQUIRE(hex.volume() == Approx(tets24Volume).epsilon(epsilon));
+    REQUIRE(hex.volume() ==
+            Approx(wedges[0].volume() + wedges[1].volume()).epsilon(epsilon));
   }
-}
+}