api_device_grid_search.hpp | Ansys Developer Portal

#pragma once
 
#include <rocky20/core/utils/rocky_math.hpp>
 
#include <cstddef>
#include <vector_types.h>
 
using namespace rocky::vector;
 
template <typename NeighborhoodOperator>
__host__ __device__ size_t iterate_over_grid_cells(
    NeighborhoodOperator& neighborhood,
    const size_t* grid_starting_index,
    const size_t* grid_to_original_map,
    const longlong3& block_start,
    const longlong3& block_end,
    const longlong3& grid_size,
    const double3& home_position,
    const size_t home_entity_id,
    const double search_radius,
    const bool not_in_block)
{
    const double squared_search_radius = search_radius * search_radius;
    const size_t grid_volume = grid_size.x * grid_size.y * grid_size.z;
    const size_t near_count = neighborhood.get_near_count();
    size_t n_thread_vicinities = 0;
    for (size_t x = block_start.x; x < block_end.x; ++x)
    {
        const size_t x_offset = x * grid_size.y * grid_size.z;
        for (size_t y = block_start.y; y < block_end.y; ++y)
        {
            const size_t xy_offset = x_offset + (y * grid_size.z);
            const size_t first_cell_id = xy_offset + block_start.z;
            const size_t last_cell_id = xy_offset + block_end.z;
            const size_t index_start = grid_starting_index[first_cell_id];
            const size_t index_end =
                (last_cell_id < grid_volume) ? grid_starting_index[last_cell_id] : near_count;
 
            for (size_t prt_id_grid = index_start; prt_id_grid < index_end; ++prt_id_grid)
            {
                const size_t near_entity_id = grid_to_original_map[prt_id_grid];
                if (neighborhood.is_valid(home_entity_id, near_entity_id, not_in_block))
                {
                    const double3 near_position = neighborhood.get_near_centroid(near_entity_id);
                    const double squared_distance = rocky::vector::dot(home_position - near_position);
                    if (squared_distance <= squared_search_radius)
                    {
                        neighborhood.operate(home_entity_id, near_entity_id);
                        n_thread_vicinities++;
                    }
                }
            }
        }
    }
    return n_thread_vicinities;
}
 
template <typename NeighborhoodOperator, typename GridStartingIndex>
__host__ __device__ size_t iterate_over_grid_cells_2(
    NeighborhoodOperator& neighborhood,
    const GridStartingIndex grid_starting_index,
    const unsigned int* __restrict__ grid_to_original_map,
    const uint3& block_start,
    const uint3& block_end,
    const uint3& grid_size,
    const unsigned int home_entity_id,
    const bool not_in_block)
{
    const size_t grid_volume = size_t(grid_size.x) * grid_size.y * grid_size.z;
    const unsigned int near_count = neighborhood.get_near_count();
    unsigned int n_thread_vicinities = 0;
    for (unsigned int x = block_start.x; x < block_end.x; ++x)
    {
        const size_t x_offset = size_t(x) * grid_size.y * grid_size.z;
        for (unsigned int y = block_start.y; y < block_end.y; ++y)
        {
            const size_t xy_offset = x_offset + (size_t(y) * grid_size.z);
            const size_t first_cell_id = xy_offset + block_start.z;
            const size_t last_cell_id = xy_offset + block_end.z;
            unsigned int index_start, index_end;
            if constexpr (std::is_pointer_v<GridStartingIndex>)
            {
                index_start = grid_starting_index[first_cell_id];
                index_end =
                    (last_cell_id < grid_volume) ? grid_starting_index[last_cell_id] : near_count;
            }
            else
            {
                grid_starting_index.get_range(
                    first_cell_id, last_cell_id, index_start, index_end, grid_volume);
            }
 
            for (unsigned int prt_id_grid = index_start; prt_id_grid < index_end; ++prt_id_grid)
            {
                const unsigned int near_entity_id = grid_to_original_map[prt_id_grid];
                if (neighborhood.is_valid(home_entity_id, near_entity_id, not_in_block))
                {
                    neighborhood.operate(home_entity_id, near_entity_id);
                    n_thread_vicinities++;
                }
            }
        }
    }
    return n_thread_vicinities;
}
 
template <typename Backend, typename NeighborhoodOperator>
__host__ __device__ void travel_vicinities_outsider(
    NeighborhoodOperator& neighborhood,
    const size_t* grid_starting_index,
    const size_t* grid_to_original_map,
    const double cell_size,
    const double search_radius,
    const longlong3 grid_size,
    double3 grid_start,
    size_t* n_vicinities)
{
    auto clamp_to_grid = [&grid_size](const longlong3& position) -> longlong3 {
        constexpr longlong3 origin { 0, 0, 0 };
        return clamp_vector(position, origin, grid_size);
    };
 
    // reduce grid start by an epsilon to avoid summing it on every position calculation
    grid_start -= (cell_size * 0.01);
    const double cell_size_inverse = 1.0 / cell_size;
    size_t n_thread_vicinities = 0;
    int start, step, stop;
 
    Backend::ThreadIteration(start, step, stop, neighborhood.get_home_count());
    for (size_t home_entity_id = start; home_entity_id < stop; home_entity_id += step)
    {
        const double3 home_position = neighborhood.get_home_centroid(home_entity_id);
        const auto home_grid_position { rocky::vector::cast<longlong3>(
            (home_position - grid_start) * cell_size_inverse) };
 
        const longlong3 block_start = clamp_to_grid(home_grid_position - 2);
        const longlong3 block_end = clamp_to_grid(home_grid_position + 3);
 
        n_thread_vicinities += iterate_over_grid_cells(
            neighborhood,
            grid_starting_index,
            grid_to_original_map,
            block_start,
            block_end,
            grid_size,
            home_position,
            home_entity_id,
            search_radius,
            false);
    }
    if (n_vicinities != nullptr && n_thread_vicinities > 0)
    {
        Backend::AtomicAdd(n_vicinities, n_thread_vicinities);
    }
}
 
template <typename Backend, typename NeighborhoodOperator>
__host__ __device__ void travel_vicinities_insider(
    NeighborhoodOperator& neighborhood,
    const size_t* grid_starting_index,
    const size_t* grid_to_original_map,
    const double cell_size,
    const double search_radius,
    const longlong3 grid_size,
    double3 grid_start,
    size_t* n_vicinities)
{
    // reduce grid start by an epsilon to avoid summing it on every position calculation
    grid_start -= (cell_size * 0.01);
    const double cell_size_inverse = 1.0 / cell_size;
    size_t n_thread_vicinities = 0;
    int start, step, stop;
 
    double3 home_position;
    size_t home_entity_id;
    auto dispatch_on_block = [&](const longlong3& _block_start,
                                 const longlong3& _block_end,
                                 const bool different_cells = true) -> void {
        n_thread_vicinities += iterate_over_grid_cells(
            neighborhood,
            grid_starting_index,
            grid_to_original_map,
            _block_start,
            _block_end,
            grid_size,
            home_position,
            home_entity_id,
            search_radius,
            different_cells);
    };
 
    Backend::ThreadIteration(start, step, stop, neighborhood.get_home_count());
    for (size_t prt_id_grid = start; prt_id_grid < stop; prt_id_grid += step)
    {
        home_entity_id = grid_to_original_map[prt_id_grid];
        if (home_entity_id == INVALID_ENTITY_ID)
        {
            continue;
        }
 
        home_position = neighborhood.get_home_centroid(home_entity_id);
        const auto home_coordinate { rocky::vector::cast<longlong3>(
            (home_position - grid_start) * cell_size_inverse) };
 
        const longlong3 block_start { max_vector({ 0ll, 0ll, 0ll }, home_coordinate - 2) };
        const longlong3 block_diag { min_vector(grid_size, home_coordinate + 1) };
        const longlong3 block_end { min_vector(grid_size, home_coordinate + 3) };
 
        // xyz block above
        dispatch_on_block({ block_diag.x, block_start.y, block_start.z }, block_end);
        // right side of yz plane
        dispatch_on_block(
            { home_coordinate.x, block_diag.y, block_start.z },
            { block_diag.x, block_end.y, block_end.z });
        // all the way back to z end on this line
        dispatch_on_block(
            { home_coordinate.x, home_coordinate.y, block_diag.z },
            { block_diag.x, block_diag.y, block_end.z });
        // same grid cell. Can't skip element ID check to avoid duplicate vicinities.
        dispatch_on_block(home_coordinate, block_diag, false);
    }
    if (n_vicinities != nullptr && n_thread_vicinities > 0)
    {
        Backend::AtomicAdd(n_vicinities, n_thread_vicinities);
    }
}
 
template <typename Backend, typename NeighborhoodOperator, typename GridStartingIndex>
__host__ __device__ unsigned int travel_vicinities_insider_full(
    NeighborhoodOperator& neighborhood,
    const GridStartingIndex grid_starting_index,
    const unsigned int* __restrict__ grid_to_original_map,
    const float cell_size,
    const float search_radius,
    const uint3& grid_size,
    float3 grid_start)
{
    const float epsilon = cell_size * 0.01f;
    const float cell_size_inverse = 1.f / cell_size;
    const unsigned int steps = std::ceil(search_radius / cell_size);
    unsigned int n_thread_vicinities = 0;
    int start, step, stop;
 
    unsigned int home_entity_id;
    auto dispatch_on_block = [&](const uint3& _block_start,
                                 const uint3& _block_end,
                                 const bool different_cells = true) -> void {
        n_thread_vicinities += iterate_over_grid_cells_2(
            neighborhood,
            grid_starting_index,
            grid_to_original_map,
            _block_start,
            _block_end,
            grid_size,
            home_entity_id,
            different_cells);
    };
 
    Backend::ThreadIteration(start, step, stop, neighborhood.get_home_count());
    for (home_entity_id = start; home_entity_id < stop; home_entity_id += step)
    {
        neighborhood.start_element(home_entity_id);
 
        const auto home_coordinate = rocky::vector::cast<uint3>(
            ((neighborhood.get_home_centroid(home_entity_id) - grid_start) + epsilon)
            * cell_size_inverse);
 
        const uint3 block_start { max_vector(uint3 { steps, steps, steps }, home_coordinate) - steps };
        const uint3 block_diag { min_vector(grid_size, home_coordinate + 1) };
        const uint3 block_end { min_vector(grid_size, home_coordinate + (steps + 1)) };
 
        // xyz block above
        dispatch_on_block({ block_diag.x, block_start.y, block_start.z }, block_end);
        // right side of yz plane
        dispatch_on_block(
            { home_coordinate.x, block_diag.y, block_start.z },
            { block_diag.x, block_end.y, block_end.z });
        // all the way back to z end on this line
        dispatch_on_block(
            { home_coordinate.x, home_coordinate.y, block_diag.z },
            { block_diag.x, block_diag.y, block_end.z });
        // same grid cell. Can't skip element ID check to avoid duplicate vicinities.
        dispatch_on_block(home_coordinate, block_diag, false);
 
        neighborhood.finish_element(home_entity_id);
    }
    return n_thread_vicinities;
}
Rocky Solver SDK 2024 R2

api_device_grid_search

Connect with Ansys
