Skip to content

GitLab

Commit 961dd20b authored by ESTERIE Pierre's avatar ESTERIE Pierre
Browse files

Add new storage type and fix examples and units

parent cff3d39b
Expand all Hide whitespace changes
Inline Side-by-side
experimental/examples/count_kernel.cpp
View file @ 961dd20b
......@@ -3,6 +3,7 @@
#include "scalfmm/container/particle.hpp"
#include "scalfmm/container/particle_container.hpp"
//
#include "scalfmm/interpolation/grid_storage.hpp"
#include "scalfmm/meta/utils.hpp"
#include "scalfmm/tools/colorized.hpp"
#include "scalfmm/tools/fma_loader.hpp"
......@@ -147,7 +148,7 @@ auto run(const int& tree_height, const int& group_size, const bool readFile, std
using particle_type = scalfmm::container::particle<double, Dimension, double, number_of_physical_values, double, 1>;
using container_type = scalfmm::container::particle_container<particle_type>;
using position_type = typename particle_type::position_type;
using cell_type = scalfmm::component::cell<double, Dimension, number_of_physical_values, 1>;
using cell_type = scalfmm::component::cell<scalfmm::component::grid_storage<double, Dimension, number_of_physical_values, 1>>;
using leaf_type = scalfmm::component::leaf<particle_type>;
using box3_type = scalfmm::component::box<position_type>;
using group_tree_type = scalfmm::component::group_tree<cell_type, leaf_type, box3_type>;
......
experimental/examples/test-like-mrhs.cpp
View file @ 961dd20b
......@@ -73,15 +73,15 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
// ---------------------------------------
// scalfmm 3.0 tree tests and benchmarks.
// ---------------------------------------
using matrix_kernel_type = scalfmm::matrix_kernels::laplace::like_mrhs;
using interpolator_type = scalfmm::interpolation::uniform_interpolator<double, dimension, matrix_kernel_type>;
using particle_type = scalfmm::container::particle<double, dimension, double, inputs, double, outputs, std::size_t>;
using container_type = scalfmm::container::particle_container<particle_type>;
using position_type = typename particle_type::position_type;
using cell_type = scalfmm::component::cell<double, dimension, inputs, outputs, std::complex<double>>;
using cell_type = scalfmm::component::cell<typename interpolator_type::storage_type>;
using leaf_type = scalfmm::component::leaf<particle_type>;
using box_type = scalfmm::component::box<position_type>;
using group_tree_type = scalfmm::component::group_tree<cell_type, leaf_type, box_type>;
using matrix_kernel_type = scalfmm::matrix_kernels::laplace::like_mrhs;
using interpolator_type = scalfmm::interpolation::uniform_interpolator<double, dimension, matrix_kernel_type>;
std::cout << scalfmm::colors::green << "Creating & Inserting " << number_of_particles
<< "particles for version .0 ...\n"
<< scalfmm::colors::reset;
......
experimental/examples/test-morton-index.cpp
View file @ 961dd20b
......@@ -30,6 +30,7 @@
#include "scalfmm/utils/sort.hpp"
#include "scalfmm/utils/tensor.hpp"
#include "scalfmm/utils/timer.hpp"
#include "scalfmm/interpolation/grid_storage.hpp"
using Value_type = double;
auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
......@@ -82,7 +83,7 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
// number_of_physical_values, Value_type, 1>;
using container_type = scalfmm::container::particle_container<particle_type>;
using position_type = typename particle_type::position_type;
using cell_type = scalfmm::component::cell<Value_type, dimension, inputs, outputs, std::complex<Value_type>>;
using cell_type = scalfmm::component::cell<scalfmm::component::grid_storage<double,dimension,inputs,outputs>>;
using leaf_type = scalfmm::component::leaf<particle_type>;
using box_type = scalfmm::component::box<position_type>;
using group_tree_type = scalfmm::component::group_tree<cell_type, leaf_type, box_type>;
......
experimental/examples/test-morton-sort.cpp
View file @ 961dd20b
......@@ -12,6 +12,7 @@
#include "scalfmm/tree/leaf.hpp"
#include "scalfmm/utils/sort.hpp"
#include "scalfmm/utils/timer.hpp"
#include "scalfmm/interpolation/grid_storage.hpp"
#include <array>
#include <bits/c++config.h>
......@@ -85,7 +86,7 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
using particle_type = scalfmm::container::particle<double, dimension, double, inputs, double, outputs, std::size_t>;
using container_type = scalfmm::container::particle_container<particle_type>;
using position_type = typename particle_type::position_type;
using cell_type = scalfmm::component::cell<double, dimension, inputs, outputs, std::complex<double>>;
using cell_type = scalfmm::component::cell<scalfmm::component::grid_storage<double,dimension,inputs,outputs>>;
using leaf_type = scalfmm::component::leaf<particle_type>;
using box_type = scalfmm::component::box<position_type>;
using group_tree_type = scalfmm::component::group_tree<cell_type, leaf_type, box_type>;
......
experimental/examples/test_dimension.cpp
View file @ 961dd20b
......@@ -152,8 +152,7 @@ auto run(const std::string& input_file, const int& tree_height, const int& group
nb_outputs_near, value_type, std::size_t>;
using container_type = scalfmm::container::particle_container<particle_type>;
using position_type = typename particle_type::position_type;
using cell_type =
scalfmm::component::cell<value_type, Dimension, nb_inputs_far, nb_outputs_far, std::complex<value_type>>;
using cell_type = scalfmm::component::cell<typename interpolator_type::storage_type>;
using leaf_type = scalfmm::component::leaf<particle_type>;
using box_type = scalfmm::component::box<position_type>;
using group_tree_type = scalfmm::component::group_tree<cell_type, leaf_type, box_type>;
......@@ -168,8 +167,6 @@ auto run(const std::string& input_file, const int& tree_height, const int& group
read_data<Dimension>(input_file, container, box_center, box_width);
time.tac();
// const std::size_t number_of_particles = std::get<0>(container->size());
// const std::size_t number_of_particles = container->number_elements();
std::cout << scalfmm::colors::green << "... Done.\n" << scalfmm::colors::reset;
std::cout << scalfmm::colors::green << "Box center = " << box_center << " box width = " << box_width
<< scalfmm::colors::reset << '\n';
......@@ -250,7 +247,7 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
using interpolator_type = scalfmm::interpolation::uniform_interpolator<double, dim, matrix_kernel_type>;
using far_field_type = scalfmm::operators::far_field_operator<interpolator_type>;
run<1, scalfmm::operators::fmm_operators<near_field_type, far_field_type>>(input_file, tree_height, group_size,
run<dim, scalfmm::operators::fmm_operators<near_field_type, far_field_type>>(input_file, tree_height, group_size,
order);
break;
}
......@@ -268,7 +265,7 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
using interpolator_type = scalfmm::interpolation::uniform_interpolator<double, dim, matrix_kernel_type>;
using far_field_type = scalfmm::operators::far_field_operator<interpolator_type>;
run<2, scalfmm::operators::fmm_operators<near_field_type, far_field_type>>(input_file, tree_height, group_size,
run<dim, scalfmm::operators::fmm_operators<near_field_type, far_field_type>>(input_file, tree_height, group_size,
order);
break;
}
......
experimental/examples/test_l2p.cpp
View file @ 961dd20b
......@@ -83,7 +83,7 @@ auto run(const std::string& title, const std::string& input_file, const int &tre
using particle_type =
scalfmm::container::particle<value_type, dimension, value_type, nb_inputs, value_type, nb_outputs_leaf, std::size_t>;
using container_type = scalfmm::container::particle_container<particle_type>;
using cell_type = scalfmm::component::cell<value_type, dimension, nb_inputs, nb_outputs, value_type/*std::complex<value_type>*/>;
using cell_type = scalfmm::component::cell<typename interpolator_type::storage_type>;
using leaf_type = scalfmm::component::leaf<particle_type>;
using position_type = typename particle_type::position_type;
//
......
experimental/examples/test_laplace_kernels.cpp
View file @ 961dd20b
......@@ -85,7 +85,7 @@ auto run(const std::string& input_file, const std::string& output_file, const in
using container_type = scalfmm::container::particle_container<particle_type>;
using position_type = typename particle_type::position_type;
using cell_type =
scalfmm::component::cell<value_type, dimension, nb_inputs_far, nb_outputs_far, std::complex<value_type>>;
scalfmm::component::cell<typename interpolator_type::storage_type>;
using leaf_type = scalfmm::component::leaf<particle_type>;
using box_type = scalfmm::component::box<position_type>;
using group_tree_type = scalfmm::component::group_tree<cell_type, leaf_type, box_type>;
......
experimental/include/scalfmm/interpolation/chebyshev.hpp
View file @ 961dd20b
......@@ -8,6 +8,7 @@
#include "scalfmm/container/point.hpp"
#include "scalfmm/container/variadic_adaptor.hpp"
#include "scalfmm/interpolation/generate_circulent.hpp"
#include "scalfmm/interpolation/grid_storage.hpp"
#include "scalfmm/interpolation/interpolator.hpp"
#include "scalfmm/interpolation/mapping.hpp"
#include "scalfmm/interpolation/traits.hpp"
......@@ -69,6 +70,7 @@ namespace scalfmm::interpolation
static constexpr std::size_t km = matrix_kernel_type::km;
using base_type = interpolator<chebyshev_interpolator<value_type, dimension, matrix_kernel_type>>;
using k_tensor_type = xt::xarray<value_type>;
using storage_type = component::grid_storage<value_type, dimension, km, kn>;
using interaction_matrix_type = xt::xtensor_fixed<k_tensor_type, xt::xshape<kn, km>>;
using base_type::base_type;
using buffer_inner_type = empty;
......
experimental/include/scalfmm/interpolation/grid_storage.hpp 0 → 100644
View file @ 961dd20b
// --------------------------------
// See LICENCE file at project root
// File : interpolation/uniform.hpp
// --------------------------------
#ifndef SCALFMM_INTERPOLATION_GRID_STORAGE_HPP
#define SCALFMM_INTERPOLATION_GRID_STORAGE_HPP
#include "scalfmm/container/variadic_adaptor.hpp"
#include <xtensor/xtensor.hpp>
#include <complex>
namespace scalfmm::component
{
// This is the base class for grid storage
// It provides the tensors for the grid storage
template<typename ValueType, std::size_t Dimension, std::size_t Inputs, std::size_t Outputs>
struct grid_storage
{
static constexpr std::size_t dimension = Dimension;
static constexpr std::size_t inputs_size = Inputs;
static constexpr std::size_t outputs_size = Outputs;
using value_type = ValueType;
using multipoles_container = xt::xtensor<value_type, dimension>;
using multipole_type = container::get_variadic_adaptor_t<multipoles_container, inputs_size>;
using local_type = container::get_variadic_adaptor_t<multipoles_container, outputs_size>;
grid_storage() = default;
grid_storage(grid_storage const&) = default;
grid_storage(grid_storage&&) noexcept = default;
inline auto operator=(grid_storage const&) -> grid_storage& = default;
inline auto operator=(grid_storage&&) noexcept -> grid_storage& = default;
~grid_storage() = default;
grid_storage(std::size_t order)
{
std::vector shape(dimension, order);
meta::for_each(m_multipoles, [&shape](auto& t) { t.resize(shape); });
meta::for_each(m_multipoles, [](auto& t) { std::fill(t.begin(), t.end(), value_type(0.)); });
meta::for_each(m_local_expansions, [&shape](auto& t) { t.resize(shape); });
meta::for_each(m_local_expansions, [](auto& t) { std::fill(t.begin(), t.end(), value_type(0.)); });
}
// Accesors to multipoles and locals
[[nodiscard]] inline auto multipoles() -> multipole_type& { return m_multipoles; }
[[nodiscard]] inline auto multipoles() const -> multipole_type const& { return m_multipoles; }
[[nodiscard]] inline auto cmultipoles() const -> multipole_type const& { return m_multipoles; }
[[nodiscard]] inline auto locals() -> local_type& { return m_local_expansions; }
[[nodiscard]] inline auto locals() const -> local_type const& { return m_local_expansions; }
[[nodiscard]] inline auto clocals() const -> local_type const& { return m_local_expansions; }
private:
multipole_type m_multipoles{};
local_type m_local_expansions{};
};
} // namespace scalfmm::component
#endif // SCALFMM_INTERPOLATION_GRID_STORAGE_HPP
experimental/include/scalfmm/interpolation/uniform.hpp
View file @ 961dd20b
......@@ -5,420 +5,7 @@
#ifndef SCALFMM_INTERPOLATION_UNIFORM_HPP
#define SCALFMM_INTERPOLATION_UNIFORM_HPP
#include "scalfmm/container/point.hpp"
#include "scalfmm/container/variadic_adaptor.hpp"
#include "scalfmm/interpolation/generate_circulent.hpp"
#include "scalfmm/interpolation/interpolator.hpp"
#include "scalfmm/interpolation/mapping.hpp"
#include "scalfmm/interpolation/traits.hpp"
#include "scalfmm/matrix_kernels/mk_common.hpp"
#include "scalfmm/meta/traits.hpp"
#include "scalfmm/meta/utils.hpp"
#include "scalfmm/simd/memory.hpp"
#include "scalfmm/simd/utils.hpp"
#include "scalfmm/tools/colorized.hpp"
#include "scalfmm/utils/math.hpp"
#include "scalfmm/utils/tensor.hpp"
#include <algorithm>
#include <array>
#include <bits/c++config.h>
#include <cassert>
#include <cmath>
#include <complex>
#include <iterator>
#include <limits>
#include <memory>
#include <tuple>
#include <type_traits>
#include <utility>
#include <xsimd/memory/xsimd_load_store.hpp>
#include <xsimd/xsimd.hpp>
#include <xtensor-fftw/basic.hpp>
#include <xtensor/xarray.hpp>
#include <xtensor/xbuilder.hpp>
#include <xtensor/xcomplex.hpp>
#include <xtensor/xio.hpp>
#include <xtensor/xmanipulation.hpp>
#include <xtensor/xmath.hpp>
#include <xtensor/xnoalias.hpp>
#include <xtensor/xslice.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xtensor_forward.hpp>
#include <xtensor/xtensor_simd.hpp>
#include <xtensor/xvectorize.hpp>
#include <xtl/xclosure.hpp>
#include <xtl/xcomplex.hpp>
namespace scalfmm::interpolation
{
// uniform_interpolator is a CRTP based class
// meaning that it implements the functions needed by the interpolator class
template<typename ValueType, std::size_t Dimension, typename FarFieldMatrixKernel>
struct uniform_interpolator : public interpolator<uniform_interpolator<ValueType, Dimension, FarFieldMatrixKernel>>
{
public:
using value_type = ValueType;
static constexpr std::size_t dimension = Dimension;
using matrix_kernel_type = FarFieldMatrixKernel;
static constexpr std::size_t kn = matrix_kernel_type::kn;
static constexpr std::size_t km = matrix_kernel_type::km;
using base_type = interpolator<uniform_interpolator<value_type, dimension, matrix_kernel_type>>;
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 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;
auto operator=(uniform_interpolator const&) -> uniform_interpolator& = default;
auto operator=(uniform_interpolator&&) noexcept -> uniform_interpolator& = default;
~uniform_interpolator() = default;
[[nodiscard]] inline auto roots_impl(std::size_t order) const
{
return xt::linspace(value_type(-1.), value_type(1), order);
}
template<typename ComputationType>
[[nodiscard]] inline auto polynomials_impl(ComputationType x, std::size_t n) const
{
using value_type = ComputationType;
// assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
const value_type two_const{2.0};
simd::compare(x);
auto order = this->order();
// Specific precomputation of scale factor
// in order to minimize round-off errors
// NB: scale factor could be hardcoded (just as the roots)
value_type scale{};
const int omn{int(order) - int(n) - 1};
if(omn % 2 != 0)
{
scale = value_type(-1.); // (-1)^(n-1-(k+1)+1)=(-1)^(omn-1)
}
else
{
scale = value_type(1.);
}
scale /= xsimd::pow(two_const, int(order) - 1) * math::factorial<value_type>(int(n)) *
math::factorial<value_type>(omn);
// compute L
value_type L{1.};
const value_type coeff = value_type(int(order) - 1);
for(std::size_t m = 0; m < order; ++m)
{
if(m != n)
{
// previous version with risks of round-off error
// L *= (x-FUnifRoots<order>::roots[m])/(FUnifRoots<order>::roots[n]-FUnifRoots<order>::roots[m]);
// new version (reducing round-off)
// regular grid on [-1,1] (h simplifies, only the size of the domain and a remains i.e. 2. and -1.)
L *= ((value_type(int(order) - 1) * (x + value_type(1.))) - (two_const * value_type(int(m))));
// L *= ((coeff * (x + value_type(1.))) - (two_const * value_type(int(m))));
}
}
L *= scale;
return L;
}
template<typename ComputationType>
[[nodiscard]] inline auto derivative_impl(ComputationType x, std::size_t n) const
{
using value_type = ComputationType;
// assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
const value_type two_const{2.0};
auto order{this->order()};
simd::compare(x);
// optimized variant
value_type NdL(0.); // init numerator
value_type DdL(1.); // init denominator
value_type tmpNdL{};
auto roots(roots_impl(order));
for(unsigned int p = 0; p < order; ++p)
{
if(p != n)
{
tmpNdL = value_type(1.);
for(unsigned int m = 0; m < order; ++m)
{
if(m != n && m != p)
{
tmpNdL *= x - roots.at(m);
}
}
NdL += tmpNdL;
DdL *= roots.at(n) - roots.at(p);
} // endif
} // p
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
{
using multipoles_container = typename Cell::multipoles_container;
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
for(std::size_t m = 0; m < km; ++m)
{
mtilde.at(m).resize(shape_out);
}
// 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);
},
scale_factor);
}
}
// postprocessing multipoles
template<typename Cell>
void apply_multipoles_postprocessing_impl(Cell& current_cell,
[[maybe_unused]] buffer_type const& products) const
{
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};
const 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};
const value_type half{0.5};
const value_type quarter{0.25};
// 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 * quarter;
}
value_type half_width{current_width * half};
// 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);
// 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)
{
// 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(
[&i](auto&&... xs) {
return container::point<value_type, dimension>{std::forward<decltype(xs)>(xs)[i]...};
},
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()) || ...))
{
// 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()));