Skip to content
GitLab
Commit 6e85e889 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files

add polynomial and derivative function for Cheb

parent c850b787
Hide whitespace changes
Inline Side-by-side
experimental/examples/test_l2p.cpp
View file @ 6e85e889
......@@ -197,7 +197,6 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
box.width(0));
auto roots2_1d = cheb.roots(static_cast<std::size_t>(order));
std::cout << "roots: " << roots2_1d << std::endl;
std::exit(EXIT_FAILURE);
///
/////////////////////////////////////////////////////////////////////////////
......@@ -226,8 +225,7 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
scalfmm::component::for_each_cell_leaf(
std::begin(gtree), std::end(gtree), static_cast<std::size_t>(tree_height),
[&pos, &order, &roots_x, &roots_y, &roots_z, &global_idx, &far_uniform](auto& cell, auto& leaf) {
// func(x) = 0.5*(x*x + y*y + z*z)
// func(x) = 0.5*(x*x + y*y + z*z)
auto func = functions::poly;
auto derx_func = functions::derx_poly;
auto dery_func = functions::dery_poly;
......@@ -242,7 +240,7 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
// Move roots inside cells x = center + width/2 roots
auto center = cell.center();
auto half_width = cell.width() * 0.5;
//
// Set the collocation points inside the cell
auto pos_x = half_width * roots_x + center[0];
auto pos_y = half_width * roots_y + center[1];
auto pos_z = half_width * roots_z + center[2];
......
experimental/include/scalfmm/interpolation/chebyshev.hpp
View file @ 6e85e889
......@@ -85,93 +85,46 @@ namespace scalfmm::interpolation
const value_type coeff = value_type(3.14159265358979323846264338327950) / (value_type(order));
// xfunction xt::cos
std::cout << " cos(n) " << (xt::cos(coeff * (xt::arange(int(order - 1), int(-1), -1) + 0.5))) << std::endl;
// const auto p = xt::eval(xt::cos(coeff * (xt::arange(int(order - 1), int(-1), -1) + 0.5)));
// return xt::cos(coeff * (xt::arange(int(order - 1), int(-1), -1) + 0.5));
return xt::linspace(value_type(-1.), value_type(1), order);
}
/**
* Chebyshev polynomials of first kind \f$ T_n(x) = \cos(n \arccos(x)) \f$
*
* @param[in] n index
* @param[in] x coordinate in [-1,1]
* @return function value
*/
template<typename ComputationType>
[[nodiscard]] inline auto polynomials_impl(ComputationType x, std::size_t order, std::size_t n) const
[[nodiscard]] inline auto polynomials_impl(ComputationType x, std::size_t /*order*/, 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);
// 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;
return value_type(cos(n * acos(x)));
}
/**
* For the derivation of the Chebyshev polynomials of first kind \f$
* \frac{\mathrm{d} T_n(x)}{\mathrm{d}x} = n U_{n-1}(x) \f$ the Chebyshev
* polynomials of second kind \f$ U_{n-1}(x) = \frac{\sin(n \arccos(x))}{\sqrt{1-x^2}} \f$
* are needed.
* @param[in] x coordinate in [-1,1]
* @param[in] order not used (just for comptibility)
* @param[in] n index
* @return function value
*/
template<typename ComputationType>
[[nodiscard]] inline auto derivative_impl(ComputationType x, std::size_t order, std::size_t n) const
[[nodiscard]] inline auto derivative_impl(ComputationType x, std::size_t /*order*/, 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);
// 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;
return value_type(n * (sin(n * acos(x))) / sqrt(1 - x * x));
}
template<typename TensorViewX, typename TensorViewY>
......@@ -186,44 +139,50 @@ namespace scalfmm::interpolation
[[nodiscard]] inline auto buffer_initialization_impl() const -> buffer_type
{
std::vector<std::size_t> shape(dimension, 2 * this->order() - 1);
shape.at(dimension - 1) = this->order();
return buffer_type(buffer_shape_type{}, buffer_inner_type(shape, 0.));
struct empty
{
};
return empty{};
// std::vector<std::size_t> shape(dimension, 2 * this->order() - 1);
// shape.at(dimension - 1) = this->order();
// return buffer_type(buffer_shape_type{}, buffer_inner_type(shape, 0.));
}
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>;
auto const& multipoles = current_cell.cmultipoles();
auto& mtilde = current_cell.transformed_multipoles();
std::vector shape_nnodes(Cell::dimension, order);
std::vector shape_transformed(Cell::dimension, std::size_t(2 * order - 1));
// Resize buffers to store fft results
for(std::size_t m = 0; m < km; ++m)
{
mtilde.at(m).resize(shape_nnodes);
}
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.));
});
auto views_in_padded = meta::apply(
padded_tensors, [&order](auto& t) { return tensor::get_view<dimension>(t, xt::range(0, order)); });
views_in_padded = meta::apply(multipoles, [](auto const& r) { return r; });
meta::repeat_indexed(
[&mtilde](std::size_t m, auto const& r) { mtilde.at(m) = xt::fftw::rfftn<dimension>(r); },
padded_tensors);
// 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>;
// auto const& multipoles = current_cell.cmultipoles();
// auto& mtilde = current_cell.transformed_multipoles();
// std::vector shape_nnodes(Cell::dimension, order);
// std::vector shape_transformed(Cell::dimension, std::size_t(2 * order - 1));
// // Resize buffers to store fft results
// for(std::size_t m = 0; m < km; ++m)
// {
// mtilde.at(m).resize(shape_nnodes);
// }
// 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.));
// });
// auto views_in_padded = meta::apply(
// padded_tensors, [&order](auto& t) { return tensor::get_view<dimension>(t, xt::range(0,
// order)); });
// views_in_padded = meta::apply(multipoles, [](auto const& r) { return r; });
// meta::repeat_indexed(
// [&mtilde](std::size_t m, auto const& r) { mtilde.at(m) = xt::fftw::rfftn<dimension>(r); },
// padded_tensors);
}
template<typename Cell, typename TupleScaleFactor>
......@@ -231,32 +190,32 @@ namespace scalfmm::interpolation
TupleScaleFactor scale_factor, std::size_t order, std::size_t tree_level,
[[maybe_unused]] buffer_type& products) const
{
auto const& transformed_multipoles = source_cell.ctransformed_multipoles();
for(std::size_t m = 0; m < km; ++m)
{
meta::repeat_indexed(
[&k, &products, &transformed_multipoles, m](std::size_t n, auto const& s_f) {
products.at(n) += transformed_multipoles.at(m) * k.at(n, m) * s_f;
},
scale_factor);
}
// auto const& transformed_multipoles = source_cell.ctransformed_multipoles();
// for(std::size_t m = 0; m < km; ++m)
// {
// meta::repeat_indexed(
// [&k, &products, &transformed_multipoles, m](std::size_t n, auto const& s_f) {
// products.at(n) += transformed_multipoles.at(m) * k.at(n, m) * s_f;
// },
// scale_factor);
// }
}
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();
meta::repeat_indexed(
[&products, order](std::size_t n, auto& expansion) {
auto view_for_update = xt::eval(tensor::get_view<dimension>(
xt::fftw::irfftn<dimension>(products.at(n), true), xt::range(0, order)));
expansion += view_for_update;
},
target_expansion);
// auto order{current_cell.order()};
// auto& target_expansion = current_cell.locals();
// meta::repeat_indexed(
// [&products, order](std::size_t n, auto& expansion) {
// auto view_for_update = xt::eval(tensor::get_view<dimension>(
// xt::fftw::irfftn<dimension>(products.at(n), true), xt::range(0, order)));
// expansion += view_for_update;
// },
// target_expansion);
}
};
......
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