paving_function.hpp

Go to the documentation of this file.

// -------------------------------------------------------//
//
// SHAMROCK code for hydrodynamics
// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
//
// -------------------------------------------------------//
 
#pragma once
 
#include "shambackends/fmt_bindings/fmt_defs.hpp"
#include "shambackends/math.hpp"
#include "shambackends/sycl.hpp"
#include "shamcomm/logs.hpp"
#include "shammath/AABB.hpp"
 
namespace shammath {
 
    template<typename Tvec, class paving_func>
    inline AABB<Tvec> f_aabb(
        const paving_func &paving, const AABB<Tvec> &aabb, int i, int j, int k) {
 
        static_assert(!paving_func::can_deform_aabb, "The paving function cannot deform the AABB");
 
        Tvec min = paving.f(aabb.lower, i, j, k);
        Tvec max = paving.f(aabb.upper, i, j, k);
 
        // if the min is greater than the max, swap them
        for (size_t d = 0; d < shambase::VectorProperties<Tvec>::dimension; ++d) {
            if (min[d] > max[d]) {
                std::swap(min[d], max[d]);
            }
        }
 
        return AABB<Tvec>{min, max};
    }
 
    template<typename Tvec, class paving_func>
    inline AABB<Tvec> f_aabb_inv(
        const paving_func &paving, const AABB<Tvec> &aabb, int i, int j, int k) {
 
        static_assert(!paving_func::can_deform_aabb, "The paving function cannot deform the AABB");
 
        Tvec min = paving.f_inv(aabb.lower, i, j, k);
        Tvec max = paving.f_inv(aabb.upper, i, j, k);
 
        // if the min is greater than the max, swap them
        for (size_t d = 0; d < shambase::VectorProperties<Tvec>::dimension; ++d) {
            if (min[d] > max[d]) {
                std::swap(min[d], max[d]);
            }
        }
 
        return AABB<Tvec>{min, max};
    }
 
    template<typename Tvec>
    struct paving_function_periodic_3d {
 
        // only translation so f(AABB) is still an AABB
        inline static constexpr bool can_deform_aabb = false;
 
        Tvec box_size; 
 
        Tvec f(Tvec x, int i, int j, int k) const { return x + box_size * Tvec{i, j, k}; }
 
        Tvec f_inv(Tvec x, int i, int j, int k) const { return x - box_size * Tvec{i, j, k}; }
 
        inline AABB<Tvec> f_aabb(const AABB<Tvec> &aabb, int i, int j, int k) const {
            return shammath::f_aabb(*this, aabb, i, j, k);
        }
 
        inline AABB<Tvec> f_aabb_inv(const AABB<Tvec> &aabb, int i, int j, int k) const {
            return shammath::f_aabb_inv(*this, aabb, i, j, k);
        }
 
        inline std::vector<std::array<int, 3>> get_paving_index_intersecting(
            AABB<Tvec> aabb) const {
 
            Tvec lower = aabb.lower / box_size;
            Tvec upper = aabb.upper / box_size;
 
            lower[0] = std::floor(lower[0]);
            lower[1] = std::floor(lower[1]);
            lower[2] = std::floor(lower[2]);
 
            i32_3 low_int = {
                static_cast<i32>(lower[0]), static_cast<i32>(lower[1]), static_cast<i32>(lower[2])};
 
            std::vector<std::array<int, 3>> indices;
 
            for (int i = low_int[0]; i <= upper[0]; ++i) {
                for (int j = low_int[1]; j <= upper[1]; ++j) {
                    for (int k = low_int[2]; k <= upper[2]; ++k) {
                        indices.push_back({i, j, k});
                    }
                }
            }
 
            return indices;
        }
    };
 
    template<typename Tvec>
    struct paving_function_general_3d {
 
        // only translation so f(AABB) is still an AABB
        inline static constexpr bool can_deform_aabb = false;
 
        using Tscal = shambase::VecComponent<Tvec>;
 
        Tvec box_size;
 
        Tvec box_center;
 
        bool is_x_periodic;
        bool is_y_periodic;
        bool is_z_periodic;
 
        Tvec f(Tvec x, int i, int j, int k) const {
            Tvec off{
                (is_x_periodic) ? 0 : (x[0] - box_center[0]) * (sham::m1pown<Tscal>(i) - 1),
                (is_y_periodic) ? 0 : (x[1] - box_center[1]) * (sham::m1pown<Tscal>(j) - 1),
                (is_z_periodic) ? 0 : (x[2] - box_center[2]) * (sham::m1pown<Tscal>(k) - 1)};
            return x + box_size * Tvec{i, j, k} + off;
        }
 
        Tvec f_inv(Tvec x, int i, int j, int k) const {
            Tvec tmp = x - box_size * Tvec{i, j, k};
            Tvec off{
                (is_x_periodic) ? 0 : (tmp[0] - box_center[0]) * (sham::m1pown<Tscal>(i) - 1),
                (is_y_periodic) ? 0 : (tmp[1] - box_center[1]) * (sham::m1pown<Tscal>(j) - 1),
                (is_z_periodic) ? 0 : (tmp[2] - box_center[2]) * (sham::m1pown<Tscal>(k) - 1)};
            return tmp + off;
        }
 
        inline AABB<Tvec> f_aabb(const AABB<Tvec> &aabb, int i, int j, int k) const {
            return shammath::f_aabb(*this, aabb, i, j, k);
        }
 
        inline AABB<Tvec> f_aabb_inv(const AABB<Tvec> &aabb, int i, int j, int k) const {
            return shammath::f_aabb_inv(*this, aabb, i, j, k);
        }
 
        inline std::vector<std::array<int, 3>> get_paving_index_intersecting(
            AABB<Tvec> aabb) const {
 
            Tvec lower = aabb.lower / box_size;
            Tvec upper = aabb.upper / box_size;
 
            lower[0] = std::floor(lower[0]);
            lower[1] = std::floor(lower[1]);
            lower[2] = std::floor(lower[2]);
 
            i32_3 low_int = {
                static_cast<i32>(lower[0]), static_cast<i32>(lower[1]), static_cast<i32>(lower[2])};
 
            std::vector<std::array<int, 3>> indices;
 
            for (int i = low_int[0]; i <= upper[0]; ++i) {
                for (int j = low_int[1]; j <= upper[1]; ++j) {
                    for (int k = low_int[2]; k <= upper[2]; ++k) {
                        indices.push_back({i, j, k});
                    }
                }
            }
 
            return indices;
        }
    };
 
    template<typename Tvec>
    struct paving_function_general_3d_shear_x {
 
        // only translation so f(AABB) is still an AABB
        inline static constexpr bool can_deform_aabb = false;
 
        using Tscal = shambase::VecComponent<Tvec>;
 
        Tvec box_size;   
        Tvec box_center; 
 
        bool is_x_periodic; 
        bool is_y_periodic; 
        bool is_z_periodic; 
 
        Tscal shear_x; 
 
        Tvec f(Tvec x, int i, int j, int k) const {
            Tvec off{
                (is_x_periodic) ? 0 : (x[0] - box_center[0]) * (sham::m1pown<Tscal>(i) - 1),
                (is_y_periodic) ? 0 : (x[1] - box_center[1]) * (sham::m1pown<Tscal>(j) - 1),
                (is_z_periodic) ? 0 : (x[2] - box_center[2]) * (sham::m1pown<Tscal>(k) - 1)};
            return x + box_size * Tvec{i, j, k} + off + shear_x * Tvec{j, 0, 0};
        }
 
        Tvec f_inv(Tvec x, int i, int j, int k) const {
            Tvec tmp = x - box_size * Tvec{i, j, k} - shear_x * Tvec{j, 0, 0};
            Tvec off{
                (is_x_periodic) ? 0 : (tmp[0] - box_center[0]) * (sham::m1pown<Tscal>(i) - 1),
                (is_y_periodic) ? 0 : (tmp[1] - box_center[1]) * (sham::m1pown<Tscal>(j) - 1),
                (is_z_periodic) ? 0 : (tmp[2] - box_center[2]) * (sham::m1pown<Tscal>(k) - 1)};
            return tmp + off;
        }
 
        inline AABB<Tvec> f_aabb(const AABB<Tvec> &aabb, int i, int j, int k) const {
            return shammath::f_aabb(*this, aabb, i, j, k);
        }
 
        inline AABB<Tvec> f_aabb_inv(const AABB<Tvec> &aabb, int i, int j, int k) const {
            return shammath::f_aabb_inv(*this, aabb, i, j, k);
        }
 
        inline std::vector<std::array<int, 3>> get_paving_index_intersecting(
            AABB<Tvec> aabb) const {
 
            Tvec lower = aabb.lower / box_size;
            Tvec upper = aabb.upper / box_size;
 
            lower[0] = std::floor(lower[0]);
            lower[1] = std::floor(lower[1]);
            lower[2] = std::floor(lower[2]);
 
            lower[0] -= 1;
            upper[0] += 1;
 
            i32_3 low_int = {
                static_cast<i32>(lower[0]), static_cast<i32>(lower[1]), static_cast<i32>(lower[2])};
 
            std::vector<std::array<int, 3>> indices;
 
            for (int i = low_int[0]; i <= upper[0]; ++i) {
                for (int j = low_int[1]; j <= upper[1]; ++j) {
                    for (int k = low_int[2]; k <= upper[2]; ++k) {
 
                        Tscal shear_x_floor = std::floor(j * shear_x / box_size[0]);
                        i32 s_x             = static_cast<i32>(shear_x_floor);
 
                        i32 i_s = i - s_x;
 
                        auto aabb_mapped = f_aabb(
                            AABB<Tvec>{box_center - box_size / 2, box_center + box_size / 2},
                            i_s,
                            j,
                            k);
 
                        if (aabb_mapped.get_intersect(aabb).is_volume_not_null()) {
                            indices.push_back({i_s, j, k});
                        }
                    }
                }
            }
 
            return indices;
        }
    };
 
} // namespace shammath