diff --git a/checks/check_1d.cpp b/checks/check_1d.cpp
index ff1438e0bb9e439f488be7129548d888c27bd0fd..522a5658786e42e886183d52fa1b93c60ff93d57 100644
--- a/checks/check_1d.cpp
+++ b/checks/check_1d.cpp
@@ -17,6 +17,7 @@
#include <algorithm>
#include <array>
+#include <random>
#include <string>
#include <tuple>
#include <utility>
@@ -35,12 +36,9 @@
#include "scalfmm/tools/fma_loader.hpp"
// Tree
#include "scalfmm/interpolation/grid_storage.hpp"
-#include "scalfmm/tree/box.hpp"
-#include "scalfmm/tree/cell.hpp"
-#include "scalfmm/tree/group_tree_view.hpp"
-#include "scalfmm/tree/leaf_view.hpp"
+#include "scalfmm/tree/tree.hpp"
//
-#include "scalfmm/lists/sequential.hpp"
+#include "scalfmm/lists/lists.hpp"
// replicated container
#include "scalfmm/meta/utils.hpp"
// Fmm operators
@@ -58,7 +56,7 @@ namespace local_args
struct matrix_kernel
{
cpp_tools::cl_parser::str_vec flags = {"--kernel", "-k"};
- std::string description = "Matrix kernels: \n 0) 1/r, 2) 1/r^2, "
+ std::string description = "Matrix kernels: \n 0) 1/x, 2) 1/r^2, "
"2) val_grad( 1/r)";
using type = int;
type def = 0;
@@ -97,6 +95,28 @@ namespace local_args
flagged
};
};
+ /**
+ * @brief
+ *
+ */
+ struct input_file
+ {
+ cpp_tools::cl_parser::str_vec flags = {"--input-file", "-fin"};
+ std::string description = "Input filename (.fma or .bfma).";
+ using type = std::string;
+ type def = "";
+ };
+ /**
+ * @brief
+ *
+ */
+ struct nbPart
+ {
+ cpp_tools::cl_parser::str_vec flags = {"--num-part", "-n"};
+ std::string description = "Number of particles).";
+ using type = int;
+ type def = 10000;
+ };
} // namespace local_args
template<typename Container>
@@ -127,7 +147,8 @@ auto read_data(const std::string& filename)
for(std::size_t idx = 0; idx < number_of_particles; ++idx)
{
loader.fillParticle(values_to_read.data(), nb_val_to_red_per_part);
- particle_type p;
+ particle_type p{};
+
std::size_t ii{0};
for(auto& e: p.position())
{
@@ -144,6 +165,61 @@ auto read_data(const std::string& filename)
return std::make_tuple(container, center, width);
}
+/// @brief generate N random particles in [-1,1]
+
+/// @tparam Container
+/// @param size
+/// @return
+template<typename Container>
+auto generate_particles(const int N)
+{
+ using container_type = Container;
+ using particle_type = typename Container::value_type;
+ using point_type = typename particle_type::position_type;
+ using value_type = typename particle_type::position_value_type;
+
+ static constexpr std::size_t dimension{particle_type::dimension};
+
+ container_type container(N);
+ point_type box_center{static_cast<value_type>(0.0)};
+ value_type box_width{static_cast<value_type>(2.0)};
+
+ std::cout << cpp_tools::colors::yellow << '\n';
+ std::cout << "[generate][nb-particles] : " << N << '\n';
+ std::cout << "[generate][box-width] : " << box_width << '\n';
+ std::cout << "[generate][box-center] : " << box_center << '\n';
+ std::cout << cpp_tools::colors::reset << '\n';
+
+ // random generator
+ // std::random_device rd;
+ // std::mt19937 gen(rd());
+ std::mt19937 gen(123);
+ std::uniform_real_distribution<value_type> dis(-1., 1.);
+ auto random_r = [&dis, &gen]() { return dis(gen); };
+
+ // inserting particles in the container
+ for(std::size_t idx = 0; idx < N; ++idx)
+ {
+ // particle_type p;
+ particle_type p{};
+ for(std::size_t d = 0; d < dimension; ++d)
+ {
+ p.position(d) = box_center[d] + random_r() * 0.5 * box_width;
+ }
+ for(auto& e: p.inputs())
+ {
+ e = random_r();
+ }
+ for(auto& e: p.outputs())
+ {
+ e = value_type(0.);
+ }
+ p.variables(idx);
+ container.insert_particle(idx, p);
+ }
+
+ return std::make_tuple(container, box_center, box_width);
+}
/**
* @brief
*
@@ -185,8 +261,9 @@ auto check_output(Container const& part, Tree const& tree) -> scalfmm::utils::ac
return error;
}
template<int Dimension, typename value_type, class FMM_OPERATOR_TYPE>
-auto run(const int& tree_height, const int& group_size, const std::size_t order, const std::string& input_file,
- const std::string& output_file, const bool check, const bool displayCells, const bool displayParticles) -> int
+auto run(const int& tree_height, const int& group_size, const std::size_t order, const int N,
+ const std::string& input_file, const std::string& output_file, const bool check, const bool displayCells,
+ const bool displayParticles) -> int
{
using near_matrix_kernel_type = typename FMM_OPERATOR_TYPE::near_field_type::matrix_kernel_type;
@@ -200,6 +277,7 @@ auto run(const int& tree_height, const int& group_size, const std::size_t order,
// variables original index
using particle_type = scalfmm::container::particle<value_type, Dimension, value_type, nb_inputs_near, value_type,
nb_outputs_near, std::size_t>;
+ // using container_type = std::vector<particle_type>;
using container_type = scalfmm::container::particle_container<particle_type>;
using position_type = typename particle_type::position_type;
@@ -221,9 +299,14 @@ auto run(const int& tree_height, const int& group_size, const std::size_t order,
position_type box_center{};
value_type box_width{};
container_type container{};
- std::tie(container, box_center, box_width) = read_data<container_type>(input_file);
-
- // read_data<Dimension>(input_file, container, box_center, box_width);
+ if(!input_file.empty())
+ {
+ std::tie(container, box_center, box_width) = read_data<container_type>(input_file);
+ }
+ else
+ {
+ std::tie(container, box_center, box_width) = generate_particles<container_type>(N);
+ }
//
////////////////////////////////////////////////////
// Build tree
@@ -261,7 +344,7 @@ auto run(const int& tree_height, const int& group_size, const std::size_t order,
//
const int separation_criterion = fmm_operator.near_field().separation_criterion();
const bool mutual = true;
- scalfmm::list::sequential::build_interaction_lists(tree, tree, separation_criterion, mutual);
+ scalfmm::list::sequential::build_interaction_lists(tree, tree, fmm_operator);
auto operator_to_proceed = scalfmm::algorithms::all;
scalfmm::algorithms::sequential::sequential(tree, fmm_operator, operator_to_proceed);
@@ -292,7 +375,7 @@ template<typename V, std::size_t D, typename MK>
using interpolator_alias = scalfmm::interpolation::interpolator<V, D, MK, option>;
template<int Dimension, typename value_type>
-auto run_general(const int& tree_height, const int& group_size, const std::size_t order, const int kernel,
+auto run_general(const int& tree_height, const int& group_size, const std::size_t order, const int kernel, const int N,
const std::string& input_file, const std::string& output_file, const bool check,
const bool displayCells, const bool displayParticles)
{
@@ -301,9 +384,9 @@ auto run_general(const int& tree_height, const int& group_size, const std::size_
// scalfmm::matrix_kernels::others::one_over_r2; using
// near_matrix_kernel_type = scalfmm::matrix_kernels::others::one_over_r2;
if(kernel == 0)
- { // one_over_r
- using far_matrix_kernel_type = scalfmm::matrix_kernels::laplace::one_over_r;
- using near_matrix_kernel_type = scalfmm::matrix_kernels::laplace::one_over_r;
+ { // one_over_x
+ using far_matrix_kernel_type = scalfmm::matrix_kernels::one_d ::one_over_x;
+ using near_matrix_kernel_type = scalfmm::matrix_kernels::one_d ::one_over_x;
using near_field_type = scalfmm::operators::near_field_operator<near_matrix_kernel_type>;
//
using interpolation_type = interpolator_alias<double, Dimension, far_matrix_kernel_type>;
@@ -311,8 +394,22 @@ auto run_general(const int& tree_height, const int& group_size, const std::size_
using fmm_operators_type = scalfmm::operators::fmm_operators<near_field_type, far_field_type>;
- return run<Dimension, value_type, fmm_operators_type>(tree_height, group_size, order, input_file, output_file,
- check, displayCells, displayParticles);
+ return run<Dimension, value_type, fmm_operators_type>(tree_height, group_size, order, N, input_file,
+ output_file, check, displayCells, displayParticles);
+ }
+ else if(kernel == -1)
+ { // one_over_x
+ using far_matrix_kernel_type = scalfmm::matrix_kernels::laplace ::one_over_r;
+ using near_matrix_kernel_type = scalfmm::matrix_kernels::laplace ::one_over_r;
+ using near_field_type = scalfmm::operators::near_field_operator<near_matrix_kernel_type>;
+ //
+ using interpolation_type = interpolator_alias<double, Dimension, far_matrix_kernel_type>;
+ using far_field_type = scalfmm::operators::far_field_operator<interpolation_type, false>;
+
+ using fmm_operators_type = scalfmm::operators::fmm_operators<near_field_type, far_field_type>;
+
+ return run<Dimension, value_type, fmm_operators_type>(tree_height, group_size, order, N, input_file,
+ output_file, check, displayCells, displayParticles);
}
else if(kernel == 1)
{ // 1/r^2
@@ -325,8 +422,8 @@ auto run_general(const int& tree_height, const int& group_size, const std::size_
using fmm_operators_type = scalfmm::operators::fmm_operators<near_field_type, far_field_type>;
- return run<Dimension, value_type, fmm_operators_type>(tree_height, group_size, order, input_file, output_file,
- check, displayCells, displayParticles);
+ return run<Dimension, value_type, fmm_operators_type>(tree_height, group_size, order, N, input_file,
+ output_file, check, displayCells, displayParticles);
}
else if(kernel == 2)
{ // val_grad_one_over_r
@@ -342,8 +439,8 @@ auto run_general(const int& tree_height, const int& group_size, const std::size_
using fmm_operators_type = scalfmm::operators::fmm_operators<near_field_type, far_field_type>;
- return run<Dimension, value_type, fmm_operators_type>(tree_height, group_size, order, input_file, output_file,
- check, displayCells, displayParticles);
+ return run<Dimension, value_type, fmm_operators_type>(tree_height, group_size, order, N, input_file,
+ output_file, check, displayCells, displayParticles);
}
else
{
@@ -359,8 +456,8 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
{ //
// Parameter handling
auto parser = cpp_tools::cl_parser::make_parser(
- cpp_tools::cl_parser::help{}, args::input_file(), args::output_file(), args::tree_height{}, args::block_size{},
- args::order{}, // args::thread_count{},
+ cpp_tools::cl_parser::help{}, local_args::nbPart{}, local_args::input_file(), args::output_file(),
+ args::tree_height{}, args::block_size{}, args::order{}, // args::thread_count{},
//, args::log_file{}, args::log_level{}
local_args::matrix_kernel{}, // periodicity
/*local_args::dimension{}, */ local_args::check{}, local_args::displayParticles{}, local_args::displayCells{});
@@ -373,12 +470,17 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
const int group_size{parser.get<args::block_size>()};
std::cout << cpp_tools::colors::blue << "<params> Group Size : " << group_size << cpp_tools::colors::reset << '\n';
- const std::string input_file{parser.get<args::input_file>()};
+ const std::string input_file{parser.get<local_args::input_file>()};
+ int N{0};
if(!input_file.empty())
{
std::cout << cpp_tools::colors::blue << "<params> Input file : " << input_file << cpp_tools::colors::reset
<< '\n';
}
+ else
+ {
+ N = parser.get<local_args::nbPart>();
+ }
const std::string output_file{parser.get<args::output_file>()};
if(!output_file.empty())
@@ -396,7 +498,7 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
//
using value_type = double;
- auto ret = run_general<dimension, value_type>(tree_height, group_size, order, kernel, input_file, output_file,
+ auto ret = run_general<dimension, value_type>(tree_height, group_size, order, kernel, N, input_file, output_file,
check, displayCells, displayParticles);
return ret;
}
\ No newline at end of file
diff --git a/include/scalfmm/matrix_kernels/scalar_kernels.hpp b/include/scalfmm/matrix_kernels/scalar_kernels.hpp
index 6c86d1b93854cd925d2cf66c922342ec989d9cc5..af4fb7ce39e27d34670cd6f4b9be7b74c9e69343 100644
--- a/include/scalfmm/matrix_kernels/scalar_kernels.hpp
+++ b/include/scalfmm/matrix_kernels/scalar_kernels.hpp
@@ -21,8 +21,123 @@
#include "scalfmm/meta/utils.hpp"
+namespace scalfmm::matrix_kernels::one_d
+{
+ /**
+ * @brief The one_over_x struct corresponds to the \f$1/x\f$ 1-d kernel.
+ *
+ * The kernel is defined as \f$k(x,y): R^{km} -> R^{kn}\f$ with \f$ kn = km = 1\f$
+ * \f$k(x,y) = \frac{1}{x - y }\f$
+ *
+ * - The kernel is homogeneous \f$ k(ax,ay) = 1/a k(x,y)\f$. The scale factor is \f$1/a\f$.
+ * - The kernel is not symmetric
+ *
+ * https://doi.org/10.1137/S0036142997329232
+ */
+ struct one_over_x
+ {
+ // Mandatory constants
+ static constexpr auto homogeneity_tag{homogeneity::homogenous}; // Specify the homogeneity of the kernel.
+ static constexpr auto symmetry_tag{symmetry::non_symmetric}; // Specify the symmetry of the kernel.
+ static constexpr std::size_t km{1}; // The number of inputs of the kernel.
+ static constexpr std::size_t kn{1}; // The number of outputs of the kernel.
+
+ // Mandatory types
+ template<typename ValueType>
+ using matrix_type = std::array<ValueType, kn * km>; // Matrix type that is used in the kernel.
+ template<typename ValueType>
+ using vector_type = std::array<ValueType, kn>; // Vector type that is used in the kernel.
+
+ /**
+ * @brief Returns the name of the kernel.
+ *
+ * @return A string representing the kernel's name.
+ */
+ const std::string name() const { return std::string("one_over_x (1-d)"); }
+
+ /**
+ * @brief Returns the mutual coefficient of size \f$kn\f$.
+ *
+ * This coefficient is used during the direct pass when the kernel is applied
+ * to compute interactions inside the leaf. Utilizing the kernel's symmetry
+ * reduces the computational complexity from \f$N^2\f$ to \f$N^2 / 2\f$ (where \f$N\f$ is
+ * the number of particles).
+ *
+ * @return A vector containing the mutual coefficients.
+ */
+ template<typename ValueType>
+ [[nodiscard]] inline constexpr auto mutual_coefficient() const
+ {
+ using decayed_type = typename std::decay_t<ValueType>;
+ return vector_type<decayed_type>{decayed_type(-1.)};
+ }
+
+ /**
+ * @brief Overload of the `()` operator to evaluate the kernel at points \f$x\f$ and \f$y\f$.
+ *
+ * @tparam ValueType1 Type of the elements in the first point.
+ * @tparam ValueType2 Type of the elements in the second point.
+ *
+ * @param x Multidimensional point (e.g., a 2D point).
+ * @param y Multidimensional point (e.g., a 2D point).
+ *
+ * @return The matrix \f$K(x, y)\f$ in a vector (row-major storage).
+ */
+ template<typename ValueType1, typename ValueType2>
+ [[nodiscard]] inline auto operator()(container::point<ValueType1, 1> const& x,
+ container::point<ValueType2, 1> const& y) const noexcept
+ -> std::enable_if_t<std::is_same_v<std::decay_t<ValueType1>, std::decay_t<ValueType2>>,
+ matrix_type<std::decay_t<ValueType1>>>
+ {
+ return evaluate(x, y);
+ }
+
+ /**
+ * @brief Evaluates the kernel at points \f$x\f$ and \f$y\f$.
+ *
+ * @tparam ValueType1 Type of the elements in the first point.
+ * @tparam ValueType2 Type of the elements in the second point.
+ *
+ * @param x Multidimensional point (e.g., a 2D point).
+ * @param y Multidimensional point (e.g., a 2D point).
+ *
+ * @return The matrix \f$K(x, y)\f$ in a vector (row-major storage).
+ */
+ template<typename ValueType1, typename ValueType2>
+ [[nodiscard]] inline auto evaluate(container::point<ValueType1, 1> const& x,
+ container::point<ValueType2, 1> const& y) const noexcept
+ -> std::enable_if_t<std::is_same_v<std::decay_t<ValueType1>, std::decay_t<ValueType2>>,
+ matrix_type<std::decay_t<ValueType1>>>
+ {
+ using decayed_type = typename std::decay_t<ValueType1>;
+ return matrix_type<decayed_type>{decayed_type(1.0) / (x[0] - y[0])};
+ }
+
+ /**
+ * @brief Returns the scale factor of the kernel.
+ *
+ * This method is used only if the kernel is homogeneous.
+ *
+ * @tparam ValueType Type of the cell width.
+ *
+ * @param cell_width The width of the cell.
+ *
+ * @return A vector representing the scale factor.
+ */
+ template<typename ValueType>
+ [[nodiscard]] inline auto scale_factor(ValueType cell_width) const noexcept
+ {
+ return vector_type<ValueType>{ValueType(1.) / cell_width};
+ }
+ static constexpr int separation_criterion{1}; // Criterion used to separate near and far field.
+ };
+} // namespace scalfmm::matrix_kernels::one_d
+
namespace scalfmm::matrix_kernels::others
{
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
/**
* @brief The one_over_r2 struct corresponds to the \f$1/r^2\f$ kernel.
*
diff --git a/include/scalfmm/matrix_kernels/template_matrix_kernel.hpp b/include/scalfmm/matrix_kernels/template_matrix_kernel.hpp
index 78a6d0dd41601d9535a3a4f56c6300d277e6e75b..47462d0eeb46541ef136b773a6a7eb1cfaab16d1 100644
--- a/include/scalfmm/matrix_kernels/template_matrix_kernel.hpp
+++ b/include/scalfmm/matrix_kernels/template_matrix_kernel.hpp
@@ -70,7 +70,7 @@ struct name
const std::string name() const { return std::string("name"); }
/**
- * @brief Returns the mutual coefficient of size \f$kn\f$.
+ * @brief Returns the mutual coefficient of size \f$kn\f$ such that \f$ f(x_i,x_j) = mutual_coefficient f(x_j,x_i)\f$.
*
* This coefficient is used during the direct pass when the kernel is applied
* to compute interactions inside the leaf. Utilizing the kernel's symmetry