Commit 728d9bbe authored by ESTERIE Pierre's avatar ESTERIE Pierre

Some more work to go, under construction 🚧

parent 0eb9b6e8
......@@ -4,6 +4,7 @@
#include <chrono>
#include <functional>
#include <iostream>
#include <utility>
#include <random>
#include <xtensor/xarray.hpp>
......@@ -13,7 +14,8 @@
#include <scalfmm/container/particle_container.hpp>
#include <scalfmm/interpolation/uniform.hpp>
#include <scalfmm/tools/colorized.hpp>
#include <scalfmm/tree/box.hpp>
#include <scalfmm/tree/leaf.hpp>
#include <scalfmm/core/operators.hpp>
int main()
{
......@@ -21,12 +23,15 @@ int main()
namespace colors = scalfmm::colors;
constexpr size_type particle_dim{3};
constexpr size_type potential_dim{3};
constexpr size_type order{6};
//constexpr size_type potential_dim{3};
constexpr size_type order{4};
using particle_type = scalfmm::container::particle<double, particle_dim, double, double, double>;
//using tuple_type = std::tuple<double,double,double>;
using point_type = scalfmm::container::point<double, particle_dim>;
using container_type = scalfmm::container::particle_container<particle_type>;
using multipole_type = scalfmm::container::variadic_container<double, double, double>;
using leaf_type = scalfmm::component::leaf<particle_type, point_type>;
using interpolator_type = scalfmm::interpolation::uniform_interpolator<double, particle_dim>;
......@@ -37,15 +42,15 @@ int main()
container_type particles_source{nb_particles};
container_type particles_target{nb_particles};
std::random_device rd;
std::mt19937 gen(rd());
//std::random_device rd;
std::mt19937 gen(33);
std::uniform_real_distribution<double> dis(-1.0, 1.0);
auto random_r = [&dis, &gen]() { return dis(gen); };
point_type center = {0., 0., 0.};
point_type center = {2.*width, 0., 0.};
auto make_particle = [&tree_height, &center, &width, &random_r]() {
auto make_particle = [&center, &width, &random_r]() {
point_type position = {random_r() * width, random_r() * width, random_r() * width};
position += center;
particle_type p(position, random_r(), random_r(), random_r());
......@@ -53,7 +58,7 @@ int main()
};
std::generate(particles_source.begin(), particles_source.end(), make_particle);
center = {2. * width, 0., 0.};
center = {0, 0., 0.};
std::generate(particles_target.begin(), particles_target.end(), make_particle);
auto start = std::chrono::high_resolution_clock::now();
......@@ -64,5 +69,13 @@ int main()
<< std::chrono::duration_cast<std::chrono::microseconds>(end - start).count() << "ms\n"
<< colors::reset;
// Construct leaf
leaf_type l(std::move(particles_source), center, width);
multipole_type multipoles{scalfmm::meta::pow(order,particle_dim)};
scalfmm::operators::apply_p2m(s, order, l, multipoles);
return 0;
}
......@@ -7,29 +7,36 @@
namespace scalfmm::concepts
{
template<typename R, typename Enabler> struct require_impl;
template<typename R> struct require_impl<R, void> { using type = R; };
template<typename Return, typename... Ts>
struct require_check : require_impl<Return,std::void_t<Ts...>>
{};
//template<typename R, typename Enabler>
//struct require_impl;
template<typename From, typename To>
using Convertible = std::enable_if_t<std::is_convertible_v<From,To>>;
//template<typename R>
//struct require_impl<R, void>
//{
// using type = R;
//};
template<typename T>
using Arithmetic = std::enable_if_t<std::is_arithmetic_v<T>>;
//template<typename Return, typename... Ts>
//struct require_check : require_impl<Return, std::void_t<Ts...>>
//{
//};
template<typename T>
using Integral = std::enable_if_t<std::is_integral_v<T>>;
//template<typename From, typename To>
//using Convertible = std::enable_if_t<std::is_convertible_v<From, To>>;
template<bool Condition>
using If = std::enable_if_t<Condition>;
//template<typename T>
//using Arithmetic = std::enable_if_t<std::is_arithmetic_v<T>>;
//template<typename T>
//using Integral = std::enable_if_t<std::is_integral_v<T>>;
} // end namespace scalfmm::concepts
//template<bool Condition>
//using If = std::enable_if_t<Condition>;
// Pseudo require macro
#define requires(...) -> typename ::concept::require_check<__VA_ARGS__>::type
} // end namespace scalfmm::concepts
#endif // SCALFMM_CONCEPTS_LIBRARY_HPP
//// Pseudo require macro
//#define requires(...)->typename ::concept ::require_check < __VA_ARGS__> ::type
#endif // SCALFMM_CONCEPTS_LIBRARY_HPP
......@@ -95,7 +95,7 @@ namespace scalfmm::container
/// Floating point type
using value_type = ValueType;
/// Position type, required by the FVariadicParticleContainer
using position_t = point<value_type, dimension>;
using position_type = point<value_type, dimension>;
/// Default constructor
particle() = default;
......@@ -121,7 +121,7 @@ namespace scalfmm::container
* that case the missing ones are zero constructed.
*/
template<typename... Ts, sfinae_check<correct_attribute_list_size<Ts...>()> = 0>
particle(const position_t& pos, Ts&&... ts)
particle(const position_type& pos, Ts&&... ts)
: particle(inria::make_index_sequence<dimension>(),
inria::make_index_sequence<std::tuple_size<tuple_data_t>::value - sizeof...(ts) - dimension>(),
pos, std::forward<Ts>(ts)...)
......@@ -139,18 +139,18 @@ namespace scalfmm::container
* \brief Position getter
*
* The position is stored in the #_data tuple, to extract it we need to
* recreate a position_t object. This is done by the position_impl() method.
* recreate a position_type object. This is done by the position_impl() method.
*
* \return A new position_t object
* \return A new position_type object
*/
[[nodiscard]] position_t position() const { return position_impl(inria::make_index_sequence<dimension>()); }
[[nodiscard]] position_type position() const { return position_impl(inria::make_index_sequence<dimension>()); }
/**
* \brief Position setter
*
* \parma pos The new position
*/
void set_position(const position_t& pos)
void set_position(const position_type& pos)
{
return set_position_impl(pos, inria::make_index_sequence<dimension>());
}
......@@ -205,7 +205,7 @@ namespace scalfmm::container
* Builds the particle, zero constructing missing arguuments.
*/
template<typename... Ts, std::size_t... Is, std::size_t... Js>
particle(inria::index_sequence<Is...>, inria::index_sequence<Js...>, const position_t& pos, Ts&&... ts)
particle(inria::index_sequence<Is...> is, inria::index_sequence<Js...> js, const position_type& pos, Ts&&... ts)
: tuple_data_t(pos[Is]..., std::forward<Ts>(ts)...,
typename std::tuple_element<dimension + Js + sizeof...(ts), tuple_data_t>::type(0)...)
{
......@@ -219,9 +219,9 @@ namespace scalfmm::container
* using arguments
*/
template<std::size_t... Is>
position_t position_impl(inria::index_sequence<Is...>) const
position_type position_impl(inria::index_sequence<Is...> is) const
{
return position_t(std::get<Is>(this->as_tuple())...);
return position_type(std::get<Is>(this->as_tuple())...);
}
/**
......@@ -233,7 +233,7 @@ namespace scalfmm::container
* \param pos new position
*/
template<std::size_t... Is>
void set_position_impl(const position_t& pos, inria::index_sequence<Is...>)
void set_position_impl(const position_type& pos, inria::index_sequence<Is...> is)
{
auto l = {std::get<Is>(this->as_tuple()) = pos[Is]...};
(void)l;
......@@ -248,7 +248,7 @@ namespace scalfmm::container
* \param pos new position
*/
template<std::size_t... Is>
types_tuple_t attributes_impl(inria::index_sequence<Is...>) const
types_tuple_t attributes_impl(inria::index_sequence<Is...> is) const
{
return types_tuple_t(std::get<Is + dimension>(*this)...);
}
......
......@@ -5,46 +5,41 @@
#ifndef SCALFMM_CONTAINER_PARTICLE_CONTAINER_HPP
#define SCALFMM_CONTAINER_PARTICLE_CONTAINER_HPP
#include <vector>
#include <tuple>
#include <vector>
#include <xsimd/config/xsimd_align.hpp>
#include <scalfmm/memory/aligned_allocator.hpp>
#include <scalfmm/container/variadic_container.hpp>
#include <scalfmm/memory/aligned_allocator.hpp>
#include <scalfmm/meta/traits.hpp>
namespace scalfmm::container
{
template<typename Particle>
class particle_container
: public variadic_container
< memory::aligned_allocator< XSIMD_DEFAULT_ALIGNMENT
, typename Particle::tuple_data_t
>
, typename Particle::tuple_data_t
>
{
public:
// base type : concatenating the particle tuple with the indexes needed
// in order to allocate the vector
using value_type = typename Particle::tuple_data_t;
using allocator_type = typename memory::aligned_allocator< XSIMD_DEFAULT_ALIGNMENT, value_type>;
using base_type = variadic_container< allocator_type, value_type>;
using size_type = typename base_type::size_type;
using difference_type = typename base_type::difference_type;
using reference = typename base_type::reference_tuple;
using const_reference = typename base_type::const_reference_tuple;
using pointer = typename base_type::pointer_tuple;
using const_pointer = typename base_type::const_pointer_tuple;
using iterator = typename base_type::iterator;
using const_iterator = typename base_type::const_iterator;
template<typename Particle>
class particle_container
: public variadic_container<typename Particle::tuple_data_t>
{
public:
// base type : concatenating the particle tuple with the indexes needed
// in order to allocate the vector
using value_type = typename Particle::tuple_data_t;
using allocator_type = typename memory::aligned_allocator<XSIMD_DEFAULT_ALIGNMENT, value_type>;
using base_type = variadic_container<value_type>;
using size_type = typename base_type::size_type;
using difference_type = typename base_type::difference_type;
using reference = typename base_type::reference_tuple;
using const_reference = typename base_type::const_reference_tuple;
using pointer = typename base_type::pointer_tuple;
using const_pointer = typename base_type::const_pointer_tuple;
using iterator = typename base_type::iterator;
using const_iterator = typename base_type::const_iterator;
using particle_type = Particle;
using particle_type = Particle;
//Forwardind constructors
using base_type::base_type;
};
}
// Forwardind constructors
using base_type::base_type;
};
} // namespace scalfmm::container
#endif // SCALFMM_CONTAINER_PARTICLE_CONTAINER_HPP
#endif // SCALFMM_CONTAINER_PARTICLE_CONTAINER_HPP
......@@ -9,6 +9,10 @@
#include <iterator>
#include <ostream>
#include <type_traits>
#include <initializer_list>
#include <algorithm>
#include <xsimd/xsimd.hpp>
#include <scalfmm/meta/const_functions.hpp>
#include <scalfmm/meta/traits.hpp>
......@@ -34,22 +38,41 @@ namespace scalfmm::container
/// Space dimension count
constexpr static const std::size_t dimension = Dim;
using base_type::base_type;
point(std::initializer_list<value_type> l)
{
std::copy(l.begin(), l.end(), this->begin());
}
point() = default;
point(point const&) = default;
point(point &&) noexcept = default;
point& operator=(point const&) = default;
point& operator=(point &&) = default;
point(value_type to_splat)
{
for(std::size_t i=0; i<dimension; ++i)
{
(*this)[i] = to_splat;
}
}
template<typename... Types>
point(Types&&... e)
: base_type{{std::forward<Types>(e)...}}
template<typename U>
point(point<U> const& other)
{
for(std::size_t i=0; i < dimension; ++i)
{
(*this)[i] = value_type(other[i]);
}
}
/** Addition assignment operator */
inline point& operator+=(const point& other)
{
for(int i = 0; auto a: other)
for(std::size_t i=0; i<dimension; ++i)
{
(*this)[i++] += a;
(*this)[i] += other[i];
}
return *this;
}
......@@ -60,11 +83,13 @@ namespace scalfmm::container
return lhs;
}
inline friend point& operator+(const point lhs, const point rhs)
inline friend point operator+(point lhs, point rhs)
{
point tmp;
for(int i = 0; auto a: tmp)
tmp = lhs[i++] + rhs[i++];
for(int i = 0; i<dimension; ++i)
{
tmp = lhs[i] + rhs[i];
}
return tmp;
}
......@@ -79,7 +104,7 @@ namespace scalfmm::container
}
/** Scalar addition */
inline friend point& operator+(point& lhs, const value_type val)
inline friend point& operator+(point& lhs, value_type val)
{
lhs += val;
return lhs;
......@@ -90,7 +115,7 @@ namespace scalfmm::container
{
for(int i = 0; auto a: other)
{
(*this)[i++] -= a;
(*this)[++i] -= a;
}
return *this;
......@@ -103,11 +128,13 @@ namespace scalfmm::container
return lhs;
}
inline friend point& operator-(const point lhs, const point rhs)
inline friend point operator-(point lhs, point rhs)
{
point tmp;
for(int i = 0; auto a: tmp)
tmp = lhs[i++] - rhs[i++];
for(int i = 0; i<dimension; ++i)
{
tmp = lhs[i] - rhs[i];
}
return tmp;
}
......@@ -122,31 +149,48 @@ namespace scalfmm::container
}
/** Scalar substraction */
inline friend point& operator-(point& lhs, const value_type val)
inline friend point& operator-(point& lhs, value_type val)
{
lhs -= val;
return lhs;
}
///** Data to data multiplication assignment */
inline point& operator*=(const point& other)
{
for(std::size_t i = 0; i < dimension; ++i)
{
(*this)[i] *= other[i];
}
return *this;
}
/** Data to data multiplication operator */
inline friend point operator*(point & lhs, const point& rhs)
{
lhs *= rhs;
return lhs;
}
/** Right scalar multiplication assignment operator */
inline point& operator*=(const value_type& val)
{
for(auto& a: *this)
for(std::size_t i=0; i<dimension; ++i)
{
a *= val;
(*this)[i] *= val;
}
return *this;
}
/** Right scalar multiplication operator */
inline friend point operator*(point lhs, const value_type val)
inline friend point operator*(point lhs, value_type const& val)
{
lhs *= val;
return lhs;
}
/** Left scalar multiplication operator */
inline friend point operator*(const value_type val, point rhs)
inline friend point operator*(value_type const& val, point rhs)
{
rhs *= val;
return rhs;
......@@ -180,59 +224,85 @@ namespace scalfmm::container
}
/** Right scalar division operator */
inline friend point operator/(point lhs, const value_type val)
inline friend point operator/(point lhs, value_type val)
{
lhs /= val;
return lhs;
}
/** Equality test operator */
template<class T>
inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
-> std::enable_if_t<meta::is_float<typename T::value_type>::value, bool>
{
auto lhs_it = lhs.begin(), rhs_it = rhs.begin();
for(std::size_t i = 0; i < Dim; ++i, ++lhs_it, ++rhs_it)
{
if(!meta::feq(*lhs_it, *rhs_it, T{1e-7}))
{
return false;
}
}
return true;
}
/** Equality test operator */
template<class T>
inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
-> std::enable_if_t<meta::is_double<typename T::value_type>::value, bool>
{
auto lhs_it = lhs.begin(), rhs_it = rhs.begin();
for(std::size_t i = 0; i < Dim; ++i, ++lhs_it, ++rhs_it)
{
if(!meta::feq(*lhs_it, *rhs_it, T{1e-13}))
{
return false;
}
}
return true;
}
/** Equality test operator */
template<class T>
inline friend auto operator==(const point<T, Dim>& lhs, const point<value_type, Dim>& rhs)
-> std::enable_if_t<std::is_integral<typename T::value_type>::value, bool>
{
auto lhs_it = lhs.begin(), rhs_it = rhs.begin();
for(std::size_t i = 0; i < Dim; i++, ++lhs_it, ++rhs_it)
{
if(*lhs_it != *rhs_it)
{
return false;
}
}
return true;
}
///** Equality test operator */
//template<class T>
//inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// -> std::enable_if_t<meta::is_float<T>::value, bool>
//{
// auto lhs_it = lhs.begin();
// auto rhs_it = rhs.begin();
// const T p{1e-7};
// for(std::size_t i = 0; i < Dim; ++i, ++lhs_it, ++rhs_it)
// {
// if(!meta::feq(*lhs_it, *rhs_it, p))
// {
// return false;
// }
// }
// return true;
//}
///** Equality test operator */
//template<class T>
//inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// -> std::enable_if_t<meta::is_double<T>::value, bool>
//{
// auto lhs_it = lhs.begin();
// auto rhs_it = rhs.begin();
// const T p{1e-13};
// for(std::size_t i = 0; i < Dim; ++i, ++lhs_it, ++rhs_it)
// {
// if(!meta::feq(*lhs_it, *rhs_it, p))
// {
// return false;
// }
// }
// return true;
//}
///** Equality test operator */
//template<class T>
//inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// -> std::enable_if_t<std::is_integral<T>::value, bool>
//{
// auto lhs_it = lhs.begin();
// auto rhs_it = rhs.begin();
// for(std::size_t i = 0; i < Dim; i++, ++lhs_it, ++rhs_it)
// {
// if(*lhs_it != *rhs_it)
// {
// return false;
// }
// }
// return true;
//}
///** Equality test operator */
//template<class T>
//inline friend auto operator==(const point<T, Dim>& lhs, const point<T, Dim>& rhs)
// -> std::enable_if_t<meta::is_simd<T>::value, bool>
//{
// auto lhs_it = lhs.begin();
// auto rhs_it = rhs.begin();
// for(std::size_t i = 0; i < Dim; i++, ++lhs_it, ++rhs_it)
// {
// auto different = *lhs_it != *rhs_it;
// if(!xsimd::any(different))
// {
// return false;
// }
// }
// return true;
//}
/** Non equality test operator */
inline friend auto operator!=(const point& lhs, const point& rhs) { return !(lhs == rhs); }
......@@ -242,7 +312,9 @@ namespace scalfmm::container
{
os << "[";
for(auto it = pos.begin(); it != pos.end() - 1; it++)
{
os << *it << ", ";
}
os << pos.back() << "]";
return os;
}
......
This source diff could not be displayed because it is too large. You can view the blob instead.
// --------------------------------
// See LICENCE file at project root
// File : core/operators.hpp
// --------------------------------
#ifndef SCALFMM_CORE_OPERATORS_HPP
#define SCALFMM_CORE_OPERATORS_HPP
#include <tuple>
#include <utility>
#include <xsimd/xsimd.hpp>
#include <xtensor/xarray.hpp>
#include <scalfmm/container/point.hpp>
#include <scalfmm/interpolation/mapping.hpp>
#include <scalfmm/meta/utils.hpp>
namespace scalfmm::simd
{
namespace impl
{
struct noop_f
{
template<typename... F>
noop_f(F... fs)
{
}
};
template<typename Position, typename TuplePtr, std::size_t... Is>
constexpr inline auto load_position(TuplePtr const& particle_ptrs, std::index_sequence<Is...> s)
{
return Position{xsimd::load_aligned(std::get<Is>(particle_ptrs))...};
}
template<typename Position, typename TuplePtr, std::size_t... Is>
constexpr inline void store_position(TuplePtr const& particle_ptrs, Position const& src,
std::index_sequence<Is...> s)
{
noop_f{(xsimd::store_aligned(std::get<Is>(particle_ptrs), std::get<Is>(src)), 0)...};
}
template<typename F, typename Container, typename... Types, std::size_t... Is>
constexpr inline auto apply_f(std::index_sequence<Is...> s, F&& f, Container&& input, Types&&... ts)
{
using return_type = typename std::remove_reference<Container>::type;
return return_type{{f(std::get<Is>(std::forward<Container>(input)), std::forward<Types>(ts)...)...}};
}
template<typename Container, std::size_t... ISs>
constexpr inline auto generate_s(std::index_sequence<ISs...> s, Container const& cont,
std::array<std::size_t, sizeof...(ISs)> const& index)
{
return meta::multiply(std::get<ISs>(cont[index[ISs]])...);
}
} // namespace impl
template<typename Position, typename TuplePtr>
constexpr inline auto load_position(TuplePtr const& particle_ptrs)
{
return impl::load_position<Position>(particle_ptrs, std::make_index_sequence<Position::dimension>{});
}
template<typename Particle, std::size_t Index, typename TuplePtr>
constexpr inline auto load_attribute(TuplePtr const& particle_ptrs)
{
return xsimd::load_aligned(std::get<Particle::dimension + Index>(particle_ptrs));
}
template<typename Particle, typename Position, typename TuplePtr>
constexpr inline void store_position(TuplePtr& particle_ptrs, Position const& src)
{
impl::store_position(particle_ptrs, src, std::make_index_sequence<Particle::dimension>{});
}
template<typename Particle, std::size_t Index, typename TuplePtr, typename Attribute>
constexpr inline void store_attribute(TuplePtr& particle_ptrs, Attribute const& src)
{
xsimd::store_aligned(std::get<Particle::dimension>(particle_ptrs), src);
}
template<std::size_t Size, typename F, typename Container, typename... Types>
constexpr inline auto apply_f(F&& f, Container&& input, Types&&... ts)
{
return impl::apply_f(std::make_index_sequence<Size>{}, std::forward<F>(f), std::forward<Container>(input),
std::forward<Types>(ts)...);
}
template<std::size_t N, typename Container>
constexpr inline auto generate_s(Container const& cont, std::array<std::size_t, N> const& index)
{
return impl::generate_s(std::make_index_sequence<N>{}, cont, index);
}
} // namespace scalfmm::simd
namespace scalfmm::operators
{
template<typename Interpolator, typename Leaf, typename Multipoles>
inline void apply_p2m(Interpolator const& interp, std::size_t order, Leaf const& leaf, Multipoles& expansion)
{
using leaf_type = Leaf;
using particle_type = typename leaf_type::particle_type;