algo: add the possibility to evaluate smooth kernels at the same position

- update the 'p2p_inner' and 'full_direct' functions to take diagonal values into account for smooth kernels.
- add trait + constexpr member to detect smooth kernels at compile-time.
- add a new check source file to test the new features

checks/CMakeLists.txt

@@ -41,6 +41,8 @@ set(source_check_files
     check_reset.cpp
     check_points.cpp
     check_block.cpp
+
+    check_smooth_kernel.cpp
 )
 
 if(${CMAKE_PROJECT_NAME}_BUILD_PBC)

checks/check_smooth_kernel.cpp

+// scalfmm
+#include "scalfmm/algorithms/fmm.hpp"
+#include "scalfmm/algorithms/full_direct.hpp"
+#include "scalfmm/container/particle.hpp"
+#include "scalfmm/interpolation/interpolation.hpp"
+#include "scalfmm/matrix_kernels/gaussian.hpp"
+#include "scalfmm/matrix_kernels/laplace.hpp"
+#include "scalfmm/meta/utils.hpp"
+#include "scalfmm/operators/fmm_operators.hpp"
+#include "scalfmm/tree/box.hpp"
+#include "scalfmm/tree/cell.hpp"
+#include "scalfmm/tree/for_each.hpp"
+#include "scalfmm/tree/group_tree_view.hpp"
+#include "scalfmm/tree/leaf_view.hpp"
+#include "scalfmm/utils/accurater.hpp"
+
+// cpp tools
+#include <cpp_tools/cl_parser/help_descriptor.hpp>
+#include <cpp_tools/cl_parser/tcli.hpp>
+
+// STL
+#include <random>
+#include <vector>
+
+namespace local_args
+{
+    struct tree_height
+    {
+        cpp_tools::cl_parser::str_vec flags = {"--tree-height", "-th"};
+        std::string description = "Height of the tree.";
+        using type = std::size_t;
+        std::string input_hint = "int"; /*!< The input hint */
+        type def = 4;
+    };
+
+    struct nb_particles
+    {
+        cpp_tools::cl_parser::str_vec flags = {"--N", "--number-particles"};
+        std::string description = "Numbre of particles to generate";
+        using type = std::size_t;
+        std::string input_hint = "int"; /*!< The input hint */
+        type def = 1000;
+    };
+
+    struct order
+    {
+        cpp_tools::cl_parser::str_vec flags = {"--order", "-o"};
+        std::string description = "Order of the approximation.";
+        using type = std::size_t;
+        std::string input_hint = "int"; /*!< The input hint */
+        type def = 4;
+    };
+
+    template<typename T>
+    std::ostream& operator<<(std::ostream& os, const std::vector<T>& vec)
+    {
+        os << "{ ";
+        for(auto el: vec)
+        {
+            os << el << ' ';
+        }
+        os << "}";
+        return os;
+    }
+
+    template<typename ParserType, typename T, typename... Args>
+    auto print_parameters(ParserType const& parser, T&& value, Args&&... args) -> void
+    {
+        T dummy;
+        std::cout << cpp_tools::colors::cyan;
+        std::cout << std::boolalpha;
+        std::cout << "[param] " << std::left << std::setw(24) << dummy.flags[0] << " = " << parser.template get<T>()
+                  << "\n";
+        if constexpr(sizeof...(args) > 0)
+        {
+            print_parameters(parser, std::forward<Args>(args)...);
+        }
+        std::cout << cpp_tools::colors::reset;
+    }
+}   // namespace local_args
+
+template<typename ValueType, std::size_t Dimension, typename MatrixKernelType, typename ParserType>
+auto run(ParserType const& parser) -> void
+{
+    using value_type = ValueType;
+    static constexpr std::size_t dimension = Dimension;
+    using matrix_kernel_type = MatrixKernelType;
+
+    static constexpr std::size_t nb_inputs{matrix_kernel_type::km};
+    static constexpr std::size_t nb_outputs{matrix_kernel_type::kn};
+
+    using particle_type =
+      scalfmm::container::particle<value_type, dimension, value_type, nb_inputs, value_type, nb_outputs, std::size_t>;
+    using position_type = typename particle_type::position_type;
+    using box_type = scalfmm::component::box<position_type>;
+
+    using container_type = std::vector<particle_type>;
+
+    // interpolation types
+    using near_field_type = scalfmm::operators::near_field_operator<matrix_kernel_type>;
+    using interpolator_type = scalfmm::interpolation::interpolator<value_type, dimension, matrix_kernel_type,
+                                                                   scalfmm::options::uniform_<scalfmm::options::fft_>>;
+    using far_field_type = scalfmm::operators::far_field_operator<interpolator_type>;
+    using fmm_operator_type = scalfmm::operators::fmm_operators<near_field_type, far_field_type>;
+
+    // tree types
+    using cell_type = scalfmm::component::cell<typename interpolator_type::storage_type>;
+    using leaf_type = scalfmm::component::leaf_view<particle_type>;
+    using group_tree_type = scalfmm::component::group_tree_view<cell_type, leaf_type, box_type>;
+
+    const std::size_t nb_particles{parser.template get<local_args::nb_particles>()};
+    const std::size_t tree_height{parser.template get<local_args::tree_height>()};
+    const std::size_t order{parser.template get<local_args::order>()};
+
+    // we construct the tree
+    container_type container(nb_particles);
+
+    const value_type box_width{2.};
+    const position_type box_center(1.);
+    box_type box(box_width, box_center);
+
+    // random generator
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<value_type> dis(0.0, 2.0);
+    auto random_r = [&dis, &gen]() { return dis(gen); };
+
+    // inserting particles in the container
+    for(std::size_t idx = 0; idx < nb_particles; ++idx)
+    {
+        // particle_type p;
+        particle_type& p = container[idx];
+        for(auto& e: p.position())
+        {
+            e = random_r();
+        }
+        for(auto& e: p.inputs())
+        {
+            e = random_r();
+        }
+        for(auto& e: p.outputs())
+        {
+            e = value_type(0.);
+        }
+        p.variables(idx);
+    }
+
+    // construct the near field
+    near_field_type near_field{};
+    // a reference on the matrix_kernel of the near_field
+    auto near_mk = near_field.matrix_kernel();
+
+    // build the approximation used in the near field
+    interpolator_type interpolator(order, tree_height, box.width(0));
+    far_field_type far_field(interpolator);
+    //  construct the fmm operator
+    fmm_operator_type fmm_operator(near_field, far_field);
+
+    const std::size_t group_size{10};   // the number of cells and leaf grouped in the tree
+    group_tree_type tree(tree_height, order, box, group_size, group_size, container);
+
+    // now we have everything to call the fmm algorithm
+    scalfmm::algorithms::fmm[scalfmm::options::_s(scalfmm::options::seq)](tree, fmm_operator);
+
+    scalfmm::algorithms::full_direct(container, near_mk);
+
+    scalfmm::utils::accurater<value_type> error{};
+
+    scalfmm::component::for_each_leaf(std::cbegin(tree), std::cend(tree),
+                                      [&container, &error](auto const& leaf)
+                                      {
+                                          // loop on the particles of the leaf
+                                          for(auto const p_ref: leaf)
+                                          {
+                                              // build a particle
+                                              const auto p = typename leaf_type::const_proxy_type(p_ref);
+                                              //
+                                              const auto& idx = std::get<0>(p.variables());
+
+                                              auto const& output_ref = container[idx].outputs();
+                                              auto const& output = p.outputs();
+
+                                              for(std::size_t i{0}; i < nb_outputs; ++i)
+                                              {
+                                                  error.add(output_ref.at(i), output.at(i));
+                                              }
+                                          }
+                                      });
+
+    std::cout << cpp_tools::colors::red;
+    std::cout << error << '\n';
+    std::cout << cpp_tools::colors::reset;
+}
+
+template<typename... Args>
+auto get_parser(int argc, char* argv[], Args&&... args)
+{
+    auto parser = cpp_tools::cl_parser::make_parser(cpp_tools::cl_parser::help{}, std::forward<Args>(args)...);
+    parser.parse(argc, argv);
+    local_args::print_parameters(parser, std::forward<Args>(args)...);
+    std::cout << "\n";
+
+    return parser;
+}
+
+auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
+{
+    using value_type = double;
+    static constexpr std::size_t dimension = 2;
+
+    const std::size_t order{8};
+    const std::size_t tree_height{4};
+    const std::size_t nb_particles{10000};
+
+    // Parameter handling
+    auto parser = get_parser(argc, argv, local_args::nb_particles{}, local_args::tree_height{}, local_args::order{});
+
+    {
+        using matrix_kernel_type = scalfmm::matrix_kernels::gaussian<value_type>;
+        run<value_type, dimension, matrix_kernel_type>(parser);
+    }
+
+    {
+        using matrix_kernel_type = scalfmm::matrix_kernels::laplace::one_over_r;
+        run<value_type, dimension, matrix_kernel_type>(parser);
+    }
+
+    return 0;
+}

include/scalfmm/algorithms/full_direct.hpp

@@ -5,13 +5,13 @@
 #ifndef SCALFMM_ALGORITHMS_FULL_DIRECT_HPP
 #define SCALFMM_ALGORITHMS_FULL_DIRECT_HPP
 
+#include "scalfmm/container/iterator.hpp"
+#include "scalfmm/container/point.hpp"
+#include "scalfmm/matrix_kernels/laplace.hpp"
 #include "scalfmm/meta/traits.hpp"
 #include "scalfmm/meta/utils.hpp"
-#include <scalfmm/container/iterator.hpp>
-#include <scalfmm/container/point.hpp>
-#include <scalfmm/matrix_kernels/laplace.hpp>
-#include <scalfmm/utils/io_helpers.hpp>
-#include <scalfmm/utils/static_assert_as_exception.hpp>
+#include "scalfmm/utils/io_helpers.hpp"
+#include "scalfmm/utils/static_assert_as_exception.hpp"
 
 #include <array>
 #include <cstddef>
@@ -20,9 +20,9 @@
 #include <omp.h>
 
 using namespace scalfmm::io;
+
 namespace scalfmm::algorithms
 {
-
     /**
     * @brief compute the direct particle interactions for the built-in particle container.
     *
@@ -80,8 +80,15 @@ namespace scalfmm::algorithms
                 }
             };
 
-            compute(0, idx);
-            compute(idx + 1, particles.size());
+            if constexpr(meta::is_smooth_v<matrix_kernel_type>)
+            {
+                compute(0, particles.size());
+            }
+            else
+            {
+                compute(0, idx);
+                compute(idx + 1, particles.size());
+            }
         }
     }
 
@@ -137,8 +144,15 @@ namespace scalfmm::algorithms
                 }
             };
 
-            compute(0, idx);
-            compute(idx + 1, particles.size());
+            if constexpr(meta::is_smooth_v<matrix_kernel_type>)
+            {
+                compute(0, particles.size());
+            }
+            else
+            {
+                compute(0, idx);
+                compute(idx + 1, particles.size());
+            }
         }
     }
 
@@ -294,7 +308,6 @@ namespace scalfmm::algorithms
             compute(0, source_particles.size());
         }
     }
-
 }   // namespace scalfmm::algorithms
 
 #endif   // SCALFMM_ALGORITHMS_FULL_DIRECT_HPP

include/scalfmm/matrix_kernels/gaussian.hpp

@@ -37,6 +37,9 @@ namespace scalfmm::matrix_kernels
         template<typename ValueType>
         using vector_type = std::array<ValueType, kn>;   // Vector type that is used in the kernel.
 
+        // Optional constant
+        static constexpr bool is_smooth = true;   // Specify the kernel is smooth and can be evaluated at x == y.
+
         value_type m_coeff{value_type(1.)};   // parameter/coefficient of the kernel.
 
         /**

include/scalfmm/matrix_kernels/mk_common.hpp

@@ -35,4 +35,11 @@ namespace scalfmm::matrix_kernels
         condionally   //< conditionally convergent series
     };
 
+    /**
+     * @brief specifies if the kernel can be evaluated at x == y (compile-time)
+     */
+    struct smooth_tag
+    {
+    };
+
 }   // namespace scalfmm::matrix_kernels

include/scalfmm/meta/traits.hpp

@@ -920,6 +920,30 @@ namespace scalfmm::meta
     template<typename T>
     inline static constexpr bool is_particle_v = is_particle<T>::value;
 
+    /**
+     * @brief Primary template handles types that do not have a nested ::is_smooth member
+     *
+     * @tparam T
+     */
+    template<class, class = void>
+    struct is_smooth : std::false_type {};
+
+    /**
+     * @brief  Specialization recognizes types that do have a nested ::is_smooth member
+     *
+     * @tparam T
+     */
+    template<class T>
+    struct is_smooth<T, std::void_t<decltype(T::is_smooth)>> : std::true_type {};
+
+    /**
+     * @brief
+     *
+     * @tparam T
+     */
+    template<typename T>
+    inline static constexpr bool is_smooth_v = is_smooth<T>::value;
+
     /**
      * @brief
      *

include/scalfmm/meta/utils.hpp

@@ -152,7 +152,7 @@ namespace scalfmm::meta
 
     /**
      *  @brief   Aliases to simplify use
-     * 
+     *
      *  How to use if
      *  @code
      *     // Define a tuple of values