Fixed process of building grid

src/main_logic/main.cpp

@@ -72,7 +72,7 @@ int main()
     if(!polyhedronDispatchFailed)
     {
         while (true)
-	{
+	    {
             get_mapnodes_reflections(simulation_map, mapnodes_reflections);
 
             for (RunIterationFunc f : iteration_runners) {

src/main_logic/simulation_logic.cu

@@ -1,7 +1,5 @@
 #ifdef COMPILE_FOR_CPU
-
 #include <cmath>
-
 #endif //COMPILE_FOR_CPU
 
 #include <cstdio>

src/simulation_objects/geometric/Face.cu

@@ -1,5 +1,9 @@
 #include <utility>
 
+#ifdef COMPILE_FOR_CPU
+#include <cmath>
+#endif //COMPILE_FOR_CPU
+
 #include "Face.cuh"
 #include "../MapNode.cuh"
 #include "Polyhedron.cuh"
@@ -110,6 +114,42 @@ __host__ __device__ bool operator==(const Face &a, const Face &b)
     return true;
 }
 
+__host__ __device__ bool Face::contains_point(SpacePoint p)
+{
+    /*
+     * Calculate the area of the face as sum of areas of triangles with vertices `vertices[0]`, `vertices[i]`,
+     * `vertices[i + 1]` for all `i` > 1.
+     * WARNING: only works for convex polyhedrons.
+     */
+    double face_area_by_triangles = 0;
+    for(int i = 1; i < n_of_vertices - 1; ++i)
+    {
+        SpacePoint a = vertices[i] - vertices[0];
+        SpacePoint b = vertices[i + 1] - vertices[0];
+
+        // Magnitude of the cross product `a % b` is the area of the parallelogram, its half is the area of the triangle
+        face_area_by_triangles += get_distance(origin, a % b) / 2;
+    }
+
+    /*
+     * Calculate the area of the face as sum of areas of triangles with vertices `p` (point to be checked),
+     * `vertices[i]`, `vertices[i + 1]` for all `i`
+     */
+    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;
+
+        // Magnitude of the cross product `a % b` is the area of the parallelogram, its half is the area of the triangle
+        face_area_with_point += get_distance(origin, a % b) / 2;
+    }
+
+    // The face contains point if and only if the two areas are equal to each other
+    return std::abs(face_area_by_triangles - face_area_with_point) < eps;
+}
+
+
 __host__ __device__ bool operator!=(const Face &a, const Face &b)
 {
     return !(a == b);

src/simulation_objects/geometric/Face.cuh

@@ -33,7 +33,8 @@ public:
      *
      * @param vertices Array of polyhedron vertices that belong to the face.
      *      Must be ordered in a special way (see note below)
-     * @param n_of_vertices Number of vertices on the face
+     * @param n_of_vertices Number of elements in the `vertices` array (i.e. number of vertices, INCLUDING the repeated
+     *      one (see below))
      *
      * @note The vertices order matters: looking on a face <b>from outside</b> the polyhedron, some vertex (let's call
      * it A) must be saved to `vertices[0]`. It's neighbour clockwise - vertex B (B is next to A clockwise) must be
@@ -125,6 +126,21 @@ public:
      */
     __host__ __device__ friend bool operator==(const Face &a, const Face &b);
 
+
+    /**
+     * Checks whether the given point belongs to the face via area calculation.
+     *
+     * Checks if the given point belongs to the face by calculating actual area of the face and area of the face as if
+     * the point was on it. If the areas match (equal +- epsilon), the point is considered to belong to the face.
+     *
+     * @param p Point to be checked
+     *
+     * @returns `true` if the point belongs to the face, `false` otherwise
+     *
+     * @warning The method only works correctly if the face is a **convex** polygon
+     */
+    __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];
 }
@@ -96,14 +103,14 @@ __host__ __device__ double Polyhedron::calculate_square_of_surface()
     double square = 0;
     for(int i = 0; i < n_of_faces; ++i)
     {
-        // Cause first vertex of face repeats again in the end the condition is `j < faces[i].get_n_of_vertices() - 2`
+        // Start with the vertex `1`, not `0`, because the 0th vertex is repeated in the end of the vertices array
         for(int j = 1; j < faces[i].get_n_of_vertices() - 1; ++j)
         {
             SpacePoint a = faces[i].get_vertices()[j + 1] - faces[i].get_vertices()[0];
             SpacePoint b = faces[i].get_vertices()[j] - faces[i].get_vertices()[0];
 
             double sign_of_square = (a % b) * faces[i].get_normal();
-            sign_of_square /= abs(sign_of_square);
+            sign_of_square /= std::abs(sign_of_square);
             square += sign_of_square * (a ^ b) / 2;
         }
     }
@@ -191,6 +198,5 @@ __host__ __device__ SpacePoint get_projected_vector_end(SpacePoint vector_start,
                                                    vector_end, current_face, &intersection_edge_vertex_id);
     }
 
-    // If vector AB does not intersect any edge of face, `vector_end` equals `b`
     return vector_end;
 }