Commit 0c635542 authored by ESTERIE Pierre's avatar ESTERIE Pierre
Browse files

Add support for non_homogenous

parent 19cfb571
......@@ -52,21 +52,27 @@ namespace scalfmm::interpolation
Derived, typename std::enable_if_t<std::is_floating_point_v<typename interpolator_traits<Derived>::value_type>>>
{
public:
using derived_type = Derived;
using value_type = typename interpolator_traits<Derived>::value_type;
using size_type = std::size_t;
static constexpr std::size_t dimension = interpolator_traits<Derived>::dimension;
using far_field_matrix_kernel_type = typename interpolator_traits<Derived>::far_field_matrix_kernel_type;
static constexpr std::size_t kn = far_field_matrix_kernel_type::kn;
static constexpr std::size_t km = far_field_matrix_kernel_type::km;
static constexpr auto homogeneity_tag = far_field_matrix_kernel_type::homogeneity_tag;
static constexpr std::size_t max_number_of_cell{7};
using variadic_container_fc_type = container::get_variadic_adaptor_t<std::vector<value_type>, kn * km>;
using size_type = std::size_t;
using array_type = xt::xarray<value_type>;
using array_shape_type = typename array_type::shape_type;
template<std::size_t d>
using tensor_type = xt::xtensor<value_type, d>;
template<std::size_t d>
using tensor_shape_type = typename tensor_type<d>::shape_type;
using derived_type = Derived;
using interaction_matrix_type = typename interpolator_traits<derived_type>::interaction_matrix_type;
using k_tensor_type = typename interpolator_traits<derived_type>::k_tensor_type;
using buffer_inner_type = typename interpolator_traits<derived_type>::buffer_inner_type;
......@@ -88,9 +94,12 @@ namespace scalfmm::interpolation
, m_cell_width_extension(cell_width_extension)
, m_separation_criterion(separation_criterion)
, m_nnodes(meta::pow(order, dimension))
, m_m2l_interactions(math::pow(max_number_of_cell, dimension))
{
set_child_parent_interpolator();
generate_interactions_matrixes();
generate_interactions_matrixes(
homogeneity_tag == matrix_kernels::homogeneity::non_homogenous ? root_cell_width : value_type(2.),
tree_height);
}
interpolator(far_field_matrix_kernel_type const& far_field, size_type order, size_type tree_height = 3,
......@@ -103,9 +112,12 @@ namespace scalfmm::interpolation
, m_cell_width_extension(cell_width_extension)
, m_separation_criterion(separation_criterion)
, m_nnodes(meta::pow(order, dimension))
, m_m2l_interactions(math::pow(max_number_of_cell, dimension))
{
set_child_parent_interpolator();
generate_interactions_matrixes();
generate_interactions_matrixes(
homogeneity_tag == matrix_kernels::homogeneity::non_homogenous ? root_cell_width : value_type(2.),
tree_height);
}
[[nodiscard]] inline auto order() const noexcept -> std::size_t { return m_order; }
......@@ -122,6 +134,23 @@ namespace scalfmm::interpolation
return m_child_parent_interpolators;
}
[[nodiscard]] inline auto interactions_matrixes() noexcept -> std::vector<interaction_matrix_type>&
{
return m_interactions_matrixes;
}
[[nodiscard]] inline auto interactions_matrixes() const noexcept -> std::vector<interaction_matrix_type> const&
{
return m_interactions_matrixes;
}
[[nodiscard]] inline auto cinteractions_matrixes() const noexcept -> std::vector<interaction_matrix_type> const&
{
return m_interactions_matrixes;
}
[[nodiscard]] inline auto m2l_interactions() const noexcept { return m_m2l_interactions; }
[[nodiscard]] inline auto roots(std::size_t order) const { return this->derived_cast().roots_impl(order); }
template<typename ComputationType>
......@@ -161,17 +190,21 @@ namespace scalfmm::interpolation
auto apply_m2l(Cell const& source_cell, Cell& target_cell, std::size_t neighbor_idx, std::size_t order,
std::size_t tree_level, [[maybe_unused]] buffer_type& products) const -> void
{
if constexpr (homogeneity_tag == matrix_kernels::homogeneity::homogenous)
std::size_t level{};
if constexpr(homogeneity_tag == matrix_kernels::homogeneity::homogenous)
{
return this->derived_cast().template apply_m2l_impl<Cell>(
source_cell, target_cell, m_fc.at(neighbor_idx), m_far_field.scale_factor(target_cell.width()), order,
tree_level, products);
level = 0;
}
else
{
static_assert(homogeneity_tag != matrix_kernels::homogeneity::non_homogenous,
"M2L for non_homogenous not yet implemented");
// rool level is 0 (i.e the simulation box) first level of cell is level 1
// we start the indexing of matrixes at 0, ence the -1 on the cell level
level = tree_level-2;
}
return this->derived_cast().template apply_m2l_impl<Cell>(
source_cell, target_cell, m_interactions_matrixes.at(level * m_m2l_interactions + neighbor_idx),
m_far_field.scale_factor(target_cell.width()), order, tree_level, products);
}
template<typename Cell>
......@@ -324,80 +357,9 @@ namespace scalfmm::interpolation
}
}
template<typename TensorViewX, typename TensorViewY>
[[nodiscard]] inline auto generate_matrix_k(TensorViewX&& X, TensorViewY&& Y, std::size_t order,
far_field_matrix_kernel_type const& far_field, std::size_t n,
std::size_t m) const -> k_tensor_type
inline auto generate_interactions_matrixes(value_type width, std::size_t tree_height) -> void
{
return this->derived_cast().template generate_matrix_k_impl<TensorViewX, TensorViewY>(X, Y, order,
far_field, n, m);
}
inline auto generate_interactions_matrixes() -> void
{
// we need to keep the tensorial view
// let the same view goes down to tensorial of point
// X and Y are Td tensors storing point.
constexpr std::size_t max_number_of_cell{7};
std::size_t interactions = math::pow(max_number_of_cell, dimension); // - math::pow(3, dimension);
std::size_t nnodes = math::pow(m_order, dimension);
m_fc.resize(interactions, interaction_matrix_type(typename interaction_matrix_type::shape_type{},
k_tensor_type(std::vector(dimension, m_order)),
xt::layout_type::row_major));
// X is the [-1,1] box reference
auto X_gen = tensor::generate_meshgrid<dimension>(roots(m_order));
auto X = std::apply(
[](auto&&... xs) { return std::make_tuple(xt::eval(std::forward<decltype(xs)>(xs))...); }, X_gen);
auto X_points = xt::xarray<container::point<value_type, dimension>>::from_shape(std::get<0>(X).shape());
auto X_flatten_views = std::apply(
[](auto&&... xs) { return std::make_tuple(xt::flatten(std::forward<decltype(xs)>(xs))...); }, X);
auto X_points_flatten_views = xt::flatten(X_points);
for(std::size_t i = 0; i < nnodes; ++i)
{
X_points_flatten_views[i] = std::apply(
[&i](auto&&... xs) {
return container::point<value_type, dimension>{std::forward<decltype(xs)>(xs)[i]...};
},
X_flatten_views);
}
std::size_t flat_idx{0};
auto generate_all_interactions = [this, &flat_idx, &X_points](auto... is) {
if(((std::abs(is) > m_separation_criterion) || ...))
{
constexpr value_type unit_width{2.0};
xt::xarray<container::point<value_type, dimension>> centers(std::vector(dimension, m_order));
xt::xarray<container::point<value_type, dimension>> Y_points(std::vector(dimension, m_order));
centers.fill(container::point<value_type, dimension>({(value_type(is) * unit_width)...}));
Y_points = X_points + centers;
auto& nm_fc_tensor = m_fc.at(flat_idx);
for(std::size_t n = 0; n < kn; ++n)
{
for(std::size_t m = 0; m < km; ++m)
{
nm_fc_tensor.at(n, m) =
std::move(generate_matrix_k(X_points, Y_points, m_order, m_far_field, n, m));
}
}
}
++flat_idx;
};
std::array<int, dimension> starts{};
std::array<int, dimension> stops{};
starts.fill(-3);
stops.fill(4);
meta::looper_range<dimension>{}(generate_all_interactions, starts, stops);
return this->derived_cast().template generate_interactions_matrixes_impl(width, tree_height);
}
protected:
......@@ -424,8 +386,9 @@ namespace scalfmm::interpolation
const value_type m_cell_width_extension{};
const value_type m_separation_criterion{};
const size_type m_nnodes{};
const size_type m_m2l_interactions{};
array_type m_child_parent_interpolators{};
std::vector<interaction_matrix_type> m_fc{};
std::vector<interaction_matrix_type> m_interactions_matrixes{};
};
} // namespace scalfmm::interpolation
......
......@@ -18,6 +18,7 @@
#include "scalfmm/tools/colorized.hpp"
#include "scalfmm/utils/math.hpp"
#include "scalfmm/utils/tensor.hpp"
#include "scalfmm/matrix_kernels/mk_common.hpp"
#include <algorithm>
#include <array>
......@@ -67,11 +68,11 @@ namespace scalfmm::interpolation
using complex_type = std::complex<value_type>;
using k_tensor_type = xt::xarray<complex_type>;
using interaction_matrix_type = xt::xtensor_fixed<k_tensor_type, xt::xshape<kn, km>>;
using base_type::base_type;
using buffer_inner_type = xt::xarray<std::complex<value_type>>;
using buffer_shape_type = xt::xshape<kn>;
using buffer_type = xt::xtensor_fixed<buffer_inner_type, buffer_shape_type>;
using base_type::base_type;
uniform_interpolator() = default;
uniform_interpolator(uniform_interpolator const&) = default;
uniform_interpolator(uniform_interpolator&&) noexcept = default;
......@@ -169,15 +170,6 @@ namespace scalfmm::interpolation
return NdL / DdL;
}
template<typename TensorViewX, typename TensorViewY>
[[nodiscard]] inline auto generate_matrix_k_impl(TensorViewX&& X, TensorViewY&& Y, std::size_t order,
matrix_kernel_type const& far_field, std::size_t n,
std::size_t m) const -> k_tensor_type
{
build<value_type, dimension, dimension> build_c{};
auto C = build_c(std::forward<TensorViewX>(X), std::forward<TensorViewY>(Y), order, far_field, n, m);
return xt::fftw::rfftn<dimension>(C);
}
[[nodiscard]] inline auto buffer_initialization_impl() const -> buffer_type
{
......@@ -259,6 +251,106 @@ namespace scalfmm::interpolation
},
target_expansion);
}
template<typename TensorViewX, typename TensorViewY>
[[nodiscard]] inline auto generate_matrix_k(TensorViewX&& X, TensorViewY&& Y, std::size_t order,
matrix_kernel_type const& far_field, std::size_t n,
std::size_t m) const -> k_tensor_type
{
build<value_type, dimension, dimension> build_c{};
auto C = build_c(std::forward<TensorViewX>(X), std::forward<TensorViewY>(Y), order, far_field, n, m);
return xt::fftw::rfftn<dimension>(C);
}
inline auto generate_interactions_matrixes_impl(value_type width, std::size_t tree_height) -> void
{
std::size_t number_of_level{1};
std::size_t number_of_interactions{this->m2l_interactions()};
// here width is the root width
// so first level of cell is of width because we skip level 0 (the root) and (the first level)):
value_type current_width{width};
// we get the half width to scale the roots
if constexpr (base_type::homogeneity_tag == matrix_kernels::homogeneity::non_homogenous)
{
// (tree_heigh = 4 -> [0-3]) the bottom cells and leaves have the same symbolic level !
// but we remove level 0 (the root ie simulaton box) and level 1.
number_of_level = tree_height-2;
current_width = width*0.25;
}
value_type half_width{current_width*0.5};
// we need to keep the tensorial view
// let the same view goes down to tensorial of point
// X and Y are Td tensors storing point.
std::size_t nnodes = math::pow(this->order(), dimension);
auto& interactions_matrixes{this->interactions_matrixes()};
interactions_matrixes.resize(number_of_interactions*number_of_level, interaction_matrix_type(typename interaction_matrix_type::shape_type{},
k_tensor_type(std::vector(dimension, this->order())),
xt::layout_type::row_major));
for(std::size_t l{0}; l<number_of_level; ++l)
{
// X is the [-1,1] box reference
std::cout << "current_roots : " << half_width * this->roots(this->order()) << '\n';
auto X_gen = tensor::generate_meshgrid<dimension>(half_width * this->roots(this->order()));
auto X = std::apply(
[](auto&&... xs) { return std::make_tuple(xt::eval(std::forward<decltype(xs)>(xs))...); }, X_gen);
auto X_points = xt::xarray<container::point<value_type, dimension>>::from_shape(std::get<0>(X).shape());
auto X_flatten_views = std::apply(
[](auto&&... xs) { return std::make_tuple(xt::flatten(std::forward<decltype(xs)>(xs))...); }, X);
auto X_points_flatten_views = xt::flatten(X_points);
for(std::size_t i = 0; i < nnodes; ++i)
{
X_points_flatten_views[i] = std::apply(
[&i](auto&&... xs) {
return container::point<value_type, dimension>{std::forward<decltype(xs)>(xs)[i]...};
},
X_flatten_views);
}
std::size_t flat_idx{0};
auto generate_all_interactions = [this, &interactions_matrixes, &flat_idx, &X_points, current_width, number_of_interactions, l](auto... is) {
if(((std::abs(is) > this->separation_criterion()) || ...))
{
//constexpr value_type unit_width{2.0};
xt::xarray<container::point<value_type, dimension>> centers(
std::vector(dimension, this->order()));
xt::xarray<container::point<value_type, dimension>> Y_points(
std::vector(dimension, this->order()));
centers.fill(container::point<value_type, dimension>({(value_type(is) * current_width)...}));
Y_points = X_points + centers;
auto& nm_fc_tensor = interactions_matrixes.at(l*number_of_interactions + flat_idx);
for(std::size_t n = 0; n < kn; ++n)
{
for(std::size_t m = 0; m < km; ++m)
{
nm_fc_tensor.at(n, m) = std::move(
generate_matrix_k(X_points, Y_points, this->order(), this->matrix_kernel(), n, m));
}
}
}
++flat_idx;
};
std::array<int, dimension> starts{};
std::array<int, dimension> stops{};
starts.fill(-3);
stops.fill(4);
meta::looper_range<dimension>{}(generate_all_interactions, starts, stops);
current_width *= 0.5;
half_width *= 0.5;
}
}
};
template<typename ValueType, std::size_t Dimension, typename FarFieldMatrixKernel>
......
......@@ -50,7 +50,7 @@ namespace scalfmm
std::sort(std::begin(tuple_of_indexes), std::end(tuple_of_indexes),
[](auto& a, auto& b) { return std::get<0>(a) < std::get<0>(b); });
}
return std::forward<std::vector<indexing_structure>>(tuple_of_indexes);
return std::vector<indexing_structure>(tuple_of_indexes);
}
///
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment