SimulationMap.cu

#include <cmath>
#include <utility>

#include "SimulationMap.cuh"
#include "geometric/Polyhedron.cuh"
#include "geometric/Face.cuh"
#include "../jones_constants.hpp"
#include "../external/random_generator.cuh"


typedef bool(MapNode::*SetNodeMethod)(MapNode *);

typedef MapNode *(MapNode::*GetNodeMethod)() const;


namespace jc = jones_constants;


__device__ double const mapnode_dist = 2 * jc::speed;


__device__ SimulationMap::SimulationMap(Polyhedron *polyhedron) :
        polyhedron(polyhedron)
{
    /*
     * Maximum number of nodes
     * It must be greater than or equal to the real number of nodes (`n_of_nodes`) after creation node grid
     */
    int max_number_of_nodes = 2 * polyhedron->calculate_square_of_surface() / (mapnode_dist * mapnode_dist);

    bool create_new_nodes = true;  // New nodes are allowed to be created

    Face *start_face = &polyhedron->get_faces()[0];
    SpacePoint start_node_coordinates = (start_face->get_vertices()[0] + start_face->get_vertices()[1] +
                                         start_face->get_vertices()[2]) / 3;

    nodes = (MapNode *)malloc(sizeof(MapNode) * max_number_of_nodes);
    nodes[0] = MapNode(polyhedron, start_face, start_node_coordinates);
    n_of_nodes = 1;

    // Direction vector of the first node to the top neighbor
    SpacePoint direction_vector = relative_point_rotation(start_node_coordinates, start_face->get_vertices()[0],
                                                          start_face->get_normal(), M_PI * 2 * rand0to1()) -
                                  start_node_coordinates;

    /**
     * Array of direction vectors from nodes with the same index
     * as in `SimulationMap::nodes` array to their top neighbors
     */
    auto *nodes_directions = (SpacePoint *)malloc(sizeof(SpacePoint) * max_number_of_nodes);
    nodes_directions[0] = direction_vector * mapnode_dist / get_distance(direction_vector, origin);

    /**
     * Boolean array where i-th element tells whether the i-th face have nodes or not
     * `does_face_have_nodes[i]` value corresponds to `polyhedron->faces[i]` face
     */
    auto *does_face_have_nodes = (bool *)malloc(sizeof(bool) * polyhedron->get_n_of_faces());
    does_face_have_nodes[0] = true;
    for(int i = 1; i < polyhedron->get_n_of_faces(); ++i)
        does_face_have_nodes[i] = false;

    /**
     * Array of pointers to the `MapNode` member functions
     * Each of them returns the particular neighbor node
     * First array element corresponds to a top neighbor,
     * the following elements correspond to the following neighbors counterclockwise
     */
    GetNodeMethod get_node_neighbors[] = {
            &MapNode::get_top,
            &MapNode::get_left,
            &MapNode::get_bottom,
            &MapNode::get_right
    };

    /**
     * Array of pointers to the `MapNode` member functions
     * Each of them sets the link from current node to the particular neighbor node
     * First array element corresponds to a top neighbor,
     * the following elements correspond to the following neighbors counterclockwise
     */
    SetNodeMethod set_node_neighbors[] = {
            &MapNode::set_top,
            &MapNode::set_left,
            &MapNode::set_bottom,
            &MapNode::set_right
    };

    // Creating new nodes until it can be done, some nodes may have less neighbors than four
    for(int current_node_id = 0; current_node_id < n_of_nodes; ++current_node_id)
    {
        MapNode &current_node = nodes[current_node_id];
        double angle = 0;
        for(int i = 0; i < 4; ++i)
        {
            if((current_node.*get_node_neighbors[i])() == nullptr)
            {
                int neighbor_node_id = get_neighbor_node_id(current_node_id, nodes_directions, angle,
                                                            does_face_have_nodes, create_new_nodes);
                if(neighbor_node_id != -1)
                {
                    (current_node.*set_node_neighbors[i])(&nodes[neighbor_node_id]);
                }
            }
            angle += M_PI_2;
        }

        if(!create_new_nodes and current_node.get_face()->get_node() == nullptr)
        {
            current_node.get_face()->set_node(&current_node, polyhedron);
        }

        if(create_new_nodes && current_node_id == n_of_nodes - 1)
        {
            /*
             * Starting from this moment: setting all pointers to neighbors that were not set earlier
             * The node that will be set as the neighbor is the closest node to hypothetical coordinates of neighbor
             */

            // Returns to the beginning of the `SimulationMap::nodes` array
            current_node_id = -1;

            // All nodes were created
            create_new_nodes = false;
        }
    }

    free(nodes_directions);
    free(does_face_have_nodes);
}

__host__ __device__ SimulationMap &SimulationMap::operator=(SimulationMap &&other) noexcept
{
    if(this != &other)
    {
        swap(n_of_nodes, other.n_of_nodes);
        swap(nodes, other.nodes);
        swap(polyhedron, other.polyhedron);
    }

    return *this;
}

__host__ __device__ SimulationMap::SimulationMap(SimulationMap &&other) noexcept
{
    nodes = nullptr;

    *this = std::move(other);
}

__device__ SimulationMap::~SimulationMap()
{
    free(nodes);
}


__device__ long SimulationMap::find_face_index(Face *face) const
{
    long index = face - &polyhedron->get_faces()[0];
    if(0 <= index && index < polyhedron->get_n_of_faces())
        return index;
    return -1;
}


__device__ SpacePoint SimulationMap::count_neighbor_node_coordinates(int current_node_id, SpacePoint top_direction,
                                                                     double angle, bool do_projection) const
{
    MapNode &current_node = nodes[current_node_id];
    Face *current_face = current_node.get_face();
    SpacePoint neighbor_coordinates = relative_point_rotation(current_node.get_coordinates(),
                                                              current_node.get_coordinates() + top_direction,
                                                              current_face->get_normal(),
                                                              angle);
    if(do_projection)
    {
        neighbor_coordinates = get_projected_vector_end(current_node.get_coordinates(), neighbor_coordinates,
                                                        current_face, polyhedron);
    }
    return neighbor_coordinates;
}


__device__ int SimulationMap::find_index_of_nearest_node(SpacePoint dest) const
{
    int nearest_mapnode_id = 0;
    for(int neighbor = 0; neighbor < n_of_nodes; ++neighbor)
    {
        if(get_distance(nodes[neighbor].get_coordinates(), dest) <
           get_distance(nodes[nearest_mapnode_id].get_coordinates(), dest))
        {
            nearest_mapnode_id = neighbor;
        }
    }
    return nearest_mapnode_id;
}


__device__ void SimulationMap::set_direction_to_top_neighbor(int current_node_id, int neighbor_node_id,
                                                             SpacePoint *nodes_directions, double angle) const
{
    MapNode &neighbor_node = nodes[neighbor_node_id];
    MapNode &current_node = nodes[current_node_id];

    if(neighbor_node.get_face() == current_node.get_face())
    {
        nodes_directions[neighbor_node_id] = nodes_directions[current_node_id];
    }
    else
    {
        SpacePoint new_direction = neighbor_node.get_coordinates() -
                                   find_intersection_with_edge(current_node.get_coordinates(),
                                                               count_neighbor_node_coordinates(current_node_id,
                                                                                               nodes_directions[current_node_id],
                                                                                               angle, false),
                                                               current_node.get_face());
        new_direction = relative_point_rotation(neighbor_node.get_coordinates(),
                                                neighbor_node.get_coordinates() + new_direction,
                                                neighbor_node.get_face()->get_normal(),
                                                -angle) -
                        neighbor_node.get_coordinates();
        nodes_directions[neighbor_node_id] = new_direction * mapnode_dist / get_distance(new_direction, origin);
    }
}


__device__ int SimulationMap::get_neighbor_node_id(int current_node_id, SpacePoint *nodes_directions, double angle,
                                                   bool *does_face_have_nodes, bool create_new_nodes)
{
    Face *current_face = nodes[current_node_id].get_face();

    // Hypothetical coordinates of neighbor node counted using direction to the top neighbor and `angle`
    SpacePoint neighbor_coordinates = count_neighbor_node_coordinates(current_node_id,
                                                                      nodes_directions[current_node_id], angle,
                                                                      true);
    Face *next_face = polyhedron->find_face_by_point(neighbor_coordinates);
    int nearest_node_id = find_index_of_nearest_node(neighbor_coordinates);
    if(!create_new_nodes || (next_face == nodes[nearest_node_id].get_face() &&
                             get_distance(nodes[nearest_node_id].get_coordinates(), neighbor_coordinates) < eps))
    {
        // Neighbor node has already existed
        return nearest_node_id;
    }
    else if(current_face == next_face || !does_face_have_nodes[find_face_index(next_face)])
    {
        // Neighbor node does not exist, but it can be created
        nodes[n_of_nodes].detach_particle();
        nodes[n_of_nodes] = MapNode(polyhedron, next_face, neighbor_coordinates);

        set_direction_to_top_neighbor(current_node_id, n_of_nodes, nodes_directions, angle);

        does_face_have_nodes[find_face_index(next_face)] = true;

        n_of_nodes++;
        return n_of_nodes - 1;
    }
    return -1;
}


__device__ int SimulationMap::get_n_of_nodes() const
{
    return this->n_of_nodes;
}

__global__ void get_n_of_nodes(const SimulationMap *const simulation_map, int *return_value)
{
    STOP_ALL_THREADS_EXCEPT_FIRST;

    *return_value = simulation_map->get_n_of_nodes();
}