utils: Test for triangle normal calculation.

src/utils/include/utils/math/triangle_functions.hpp

@@ -24,7 +24,14 @@
 #include <cmath>
 
 namespace Utils {
-/** Computes the normal vector to the plane given by points P1P2P3 */
+/**
+ * @brief Computes the normal vector of a triangle.
+ *
+ * The sign convention is such that P1P2, P1P3 and
+ * the normal form a right-handed system.
+ * The normal vector is not normalized, e.g. its length
+ * is arbitrary.
+ */
 inline Vector3d get_n_triangle(const Vector3d &P1, const Vector3d &P2,
                                const Vector3d &P3) {
   auto const u = P2 - P1;

src/utils/tests/CMakeLists.txt

@@ -41,6 +41,7 @@ unit_test(NAME mask_test SRC mask_test.cpp DEPENDS utils)
 unit_test(NAME type_traits_test SRC type_traits_test.cpp DEPENDS utils)
 unit_test(NAME uniform_test SRC uniform_test.cpp DEPENDS utils)
 unit_test(NAME memcpy_archive_test SRC memcpy_archive_test.cpp DEPENDS utils)
+unit_test(NAME triangle_functions_test SRC triangle_functions_test.cpp DEPENDS utils)
 
 unit_test(NAME gather_buffer_test SRC gather_buffer_test.cpp DEPENDS utils Boost::mpi MPI::MPI_CXX NUM_PROC 4)
 unit_test(NAME scatter_buffer_test SRC scatter_buffer_test.cpp DEPENDS utils Boost::mpi MPI::MPI_CXX NUM_PROC 4)

src/utils/tests/triangle_functions_test.cpp

@@ -0,0 +1,68 @@
+/*
+ * Copyright (C) 2019 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/* Unit test for Utils triangle algorithms. */
+
+#define BOOST_TEST_MODULE Utils::triangle_functions
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include "utils/math/triangle_functions.hpp"
+
+auto const epsilon = std::numeric_limits<double>::epsilon();
+
+BOOST_AUTO_TEST_CASE(normal_) {
+  /* Non-degenerate case */
+  {
+    auto const P1 = Utils::Vector3d{0, 0, 0};
+    auto const P2 = Utils::Vector3d{0, 1, 0};
+    auto const P3 = Utils::Vector3d{0, 0, 1};
+
+    auto const normal = Utils::get_n_triangle(P1, P2, P3);
+
+    auto const P1P2 = P2 - P1;
+    auto const P1P3 = P3 - P1;
+
+    /* Normal should be orthogonal to the sides of the triangle */
+    BOOST_CHECK_SMALL(normal * P1P2, epsilon);
+    BOOST_CHECK_SMALL(normal * P1P3, epsilon);
+
+    /* Orientation */
+    BOOST_CHECK(vector_product(normal, P1P2) * P1P3 > 0.);
+    BOOST_CHECK(vector_product(P1P3, normal) * P1P2 > 0.);
+    BOOST_CHECK_SMALL(((normal + Utils::get_n_triangle(P2, P1, P3)).norm()),
+                      epsilon);
+  }
+
+  /* Degenerate case */
+  {
+    auto const P1 = Utils::Vector3d{-1, 0, 0};
+    auto const P2 = Utils::Vector3d{0, 0, 0};
+    auto const P3 = Utils::Vector3d{1, 0, 0};
+
+    auto const normal = Utils::get_n_triangle(P1, P2, P3);
+
+    auto const P1P2 = P2 - P1;
+    auto const P1P3 = P3 - P1;
+
+    /* Normal should be orthogonal to the sides of the triangle */
+    BOOST_CHECK_SMALL(normal * P1P2, epsilon);
+    BOOST_CHECK_SMALL(normal * P1P3, epsilon);
+  }
+}