Commit 83579aec authored by ESTERIE Pierre's avatar ESTERIE Pierre
Browse files

coments on uniform.hpp

parent 0c635542
......@@ -170,15 +170,19 @@ namespace scalfmm::interpolation
return NdL / DdL;
}
// Returns the buffers initialized for the optimized fft
[[nodiscard]] inline auto buffer_initialization_impl() const -> buffer_type
{
// shape for the output fft is [2*order-1, ..., order]
std::vector<std::size_t> shape(dimension, 2 * this->order() - 1);
shape.at(dimension - 1) = this->order();
// Here we return the tensors initialized with zeros
// we have kn (number of outputs) tensors
return buffer_type(buffer_shape_type{}, buffer_inner_type(shape, 0.));
}
// Preprocessing function for applying ffts on multipoles
// This function is called as soon as the multipoles are computed and available
template<typename Cell>
void apply_multipoles_preprocessing_impl(Cell& current_cell, std::size_t order) const
{
......@@ -186,23 +190,29 @@ namespace scalfmm::interpolation
using fft_type = container::get_variadic_adaptor_t<xt::xarray<typename multipoles_container::value_type>,
Cell::inputs_size>;
// Get the multipoles and transformed mutlipoles tensors
auto const& multipoles = current_cell.cmultipoles();
auto& mtilde = current_cell.transformed_multipoles();
// Input fft tensor shape i.e the multipoles padded with zeros
std::vector shape_transformed(dimension, std::size_t(2 * order - 1));
// Output fft tensor shape
std::vector<std::size_t> shape_out(shape_transformed);
shape_out.at(dimension - 1) = this->order();
// Creating padded tensors for the ffts inputs
fft_type padded_tensors;
meta::for_each(padded_tensors, [&shape_transformed](auto& t) { t.resize(shape_transformed); });
meta::for_each(padded_tensors, [](auto& t) {
std::fill(t.begin(), t.end(), typename multipoles_container::value_type(0.));
});
// Get the views of padded tensors to put the multipoles inside
// padded[O, 2*order-1] -> views_in_padded[0,order]
// The views will update the padded tensor with the multipoles values.
auto views_in_padded = meta::apply(
padded_tensors, [&order](auto& t) { return tensor::get_view<dimension>(t, xt::range(0, order)); });
// Copy the multipoles inside the padded tensors via the views
views_in_padded = meta::apply(multipoles, [](auto const& r) { return r; });
// Resize buffers to store fft results
......@@ -211,23 +221,26 @@ namespace scalfmm::interpolation
mtilde.at(m).resize(shape_out);
}
//meta::repeat_indexed(
// [&mtilde, this](std::size_t m, auto const& r) { fft.execute_plan(r, mtilde.at(m)); },
// padded_tensors);
// Applying real ffts on padded tensors -> km (inputs size) ffts
meta::repeat_indexed(
[&mtilde](std::size_t m, auto const& r) { mtilde.at(m) = xt::fftw::rfftn<dimension>(r); },
padded_tensors);
}
// m2l operator impl
template<typename Cell, typename TupleScaleFactor>
void apply_m2l_impl(Cell const& source_cell, Cell& target_cell, interaction_matrix_type const& k,
TupleScaleFactor scale_factor, std::size_t order, std::size_t tree_level,
[[maybe_unused]] buffer_type& products) const
{
// get the transformed multipoles
auto const& transformed_multipoles = source_cell.ctransformed_multipoles();
// we generate km*kn products
// loop on km
for(std::size_t m = 0; m < km; ++m)
{
// meta loop on kn
meta::repeat_indexed(
[&k, &products, &transformed_multipoles, m](std::size_t n, auto const& s_f) {
tensor::product(products.at(n), transformed_multipoles.at(m), k.at(n,m), s_f);
......@@ -236,6 +249,7 @@ namespace scalfmm::interpolation
}
}
// postprocessing multipoles
template<typename Cell>
void apply_multipoles_postprocessing_impl(Cell& current_cell,
[[maybe_unused]] buffer_type const& products) const
......@@ -243,25 +257,32 @@ namespace scalfmm::interpolation
auto order{current_cell.order()};
auto& target_expansion = current_cell.locals();
// applying kn inverse real ffts
meta::repeat_indexed(
[&products, order, this](std::size_t n, auto& expansion) {
// we get the view [0,order] from the output of the ifft.
auto view_for_update = xt::eval(tensor::get_view<dimension>(
xt::fftw::irfftn<dimension>(products.at(n), true), xt::range(0, order)));
// accumulate the results in local expansions
expansion += view_for_update;
},
target_expansion);
}
//This function generates the circulent tensors for generating interaction matrixes.
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{};
//generate the circulent tensor
auto C = build_c(std::forward<TensorViewX>(X), std::forward<TensorViewY>(Y), order, far_field, n, m);
// return the fft
return xt::fftw::rfftn<dimension>(C);
}
// This function generate the vector of interaction matrixes
inline auto generate_interactions_matrixes_impl(value_type width, std::size_t tree_height) -> void
{
std::size_t number_of_level{1};
......@@ -286,26 +307,39 @@ namespace scalfmm::interpolation
std::size_t nnodes = math::pow(this->order(), dimension);
// get the ref of the vector
auto& interactions_matrixes{this->interactions_matrixes()};
// resizing and initializing the vector
// homogenous -> only one level
// non_homogenous -> we generate interaction matrices for each tree level processed by the m2l operator
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));
// loop on levels
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';
// homogenous -> X is the [-1,1] box reference
// non_homogenous -> X is the [-cell_width_at_level/2, cell_width_at_level/2]
// X is a multidimensional grid generator returning a grid for X, another one for Y,...
// X is a multidimensional grid of points scaled on the size of cell
auto X_gen = tensor::generate_meshgrid<dimension>(half_width * this->roots(this->order()));
// we evaluate the generator
auto X = std::apply(
[](auto&&... xs) { return std::make_tuple(xt::eval(std::forward<decltype(xs)>(xs))...); }, X_gen);
// get an array of points
auto X_points = xt::xarray<container::point<value_type, dimension>>::from_shape(std::get<0>(X).shape());
//we flatten the grids
auto X_flatten_views = std::apply(
[](auto&&... xs) { return std::make_tuple(xt::flatten(std::forward<decltype(xs)>(xs))...); }, X);
//we flatten the tensor of points
auto X_points_flatten_views = xt::flatten(X_points);
// here, we reconstruct the points for the grids.
// i.e p(x,y,z)[i] <- x[i], y[i], z[i]
for(std::size_t i = 0; i < nnodes; ++i)
{
X_points_flatten_views[i] = std::apply(
......@@ -315,19 +349,25 @@ namespace scalfmm::interpolation
X_flatten_views);
}
// lambda for generation the matrixes corresponding to one interaction
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};
// constructing centers to generate Y
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)...}));
// here we fill the centers with the current loop indexes
centers.fill(container::point<value_type, dimension>({(value_type(is) * current_width)...}));
// then we calculate the Y points
Y_points = X_points + centers;
// we get the ref of interaction matrices to generate
auto& nm_fc_tensor = interactions_matrixes.at(l*number_of_interactions + flat_idx);
// and we generate each matrix needed for the product (kn*km matrices)
for(std::size_t n = 0; n < kn; ++n)
{
for(std::size_t m = 0; m < km; ++m)
......@@ -340,12 +380,16 @@ namespace scalfmm::interpolation
++flat_idx;
};
//loop range [-3,4[
std::array<int, dimension> starts{};
std::array<int, dimension> stops{};
starts.fill(-3);
stops.fill(4);
// he we expand at compile time d loops of the range
// the indices of the d loops are input parameters of the lambda generate_all_interactions
meta::looper_range<dimension>{}(generate_all_interactions, starts, stops);
// we divide the widths for the next tree level
current_width *= 0.5;
half_width *= 0.5;
}
......@@ -353,6 +397,7 @@ namespace scalfmm::interpolation
};
// traits class to register types inside the interpolator generic class
template<typename ValueType, std::size_t Dimension, typename FarFieldMatrixKernel>
struct interpolator_traits<uniform_interpolator<ValueType, Dimension, FarFieldMatrixKernel>>
{
......
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