Introduced Face::contains_point to check if point belongs to face i…

…nstead of previous approach.

This might have fixed one or many of: #29, #31, #38

Co-authored-by: Nikolay Nechaev <[email protected]>

src/simulation_objects/geometric/Face.cu

@@ -1,4 +1,9 @@
 #include <utility>
+#include <cassert>
+
+#ifdef COMPILE_FOR_CPU
+#include <cmath>
+#endif //COMPILE_FOR_CPU
 
 #include "Face.cuh"
 #include "../MapNode.cuh"
@@ -110,6 +115,26 @@ __host__ __device__ bool operator==(const Face &a, const Face &b)
     return true;
 }
 
+__host__ __device__ bool Face::contains_point(SpacePoint p)
+{
+    // For squares only!
+    assert(n_of_vertices == 4);
+    double face_area = get_distance(vertices[1], vertices[0]) * get_distance(vertices[2], vertices[1]);
+
+    double face_area_with_point = 0;
+    for(int i = 0; i < n_of_vertices; ++i)
+    {
+        SpacePoint this_vertex = vertices[i], next_vertex = vertices[(i + 1) % n_of_vertices];
+        SpacePoint a = this_vertex - p, b = next_vertex - p;
+        double triangle_area = 1. / 2 * sqrt((a * a) * (b * b) - (a * b) * (a * b));
+        face_area_with_point += triangle_area;
+    }
+
+    // The face contains point if and only if the two areas are equal to each other
+    return std::abs(face_area - face_area_with_point) < eps;
+}
+
+
 __host__ __device__ bool operator!=(const Face &a, const Face &b)
 {
     return !(a == b);

src/simulation_objects/geometric/Face.cuh

@@ -125,6 +125,9 @@ public:
      */
     __host__ __device__ friend bool operator==(const Face &a, const Face &b);
 
+
+    __host__ __device__ bool contains_point(SpacePoint p);
+
 private:
     /// Pointer to some node laying on the face
     MapNode *node;

src/simulation_objects/geometric/Polyhedron.cu

@@ -72,10 +72,17 @@ __host__ __device__ Face *Polyhedron::find_face_by_point(SpacePoint point) const
     for(int i = 0; i < n_of_faces; ++i)
     {
         Face *face = &faces[i];
+
+        /*
+        // Old version of the check. Probably will be investigated later?
         SpacePoint normal = (point - face->get_vertices()[0]) % (face->get_vertices()[1] - face->get_vertices()[0]);
         normal = normal / get_distance(normal, origin);
         if(normal * face->get_normal() >= 1 - eps)
             return face;
+        */
+
+        if(face->contains_point(point))
+            return face;
     }
     return &faces[0];
 }