From eb97ceef955fa1c5e3794c74fc9894fda1ce6f21 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 2 Jun 2026 18:15:13 +0200 Subject: Doing aabb computations --- lib/core/c++/chunk.hpp | 5 ++++ lib/core/c++/hlbm.hpp | 43 +++++++++++++++++++++++++--- lib/core/c++/particle/aabb.hpp | 49 ++++++++++++++++++++++++++++++++ lib/core/c++/particle/porosity.hpp | 57 ++++++++++++++++++++++++++++++-------- 4 files changed, 139 insertions(+), 15 deletions(-) create mode 100644 lib/core/c++/particle/aabb.hpp (limited to 'lib/core') diff --git a/lib/core/c++/chunk.hpp b/lib/core/c++/chunk.hpp index 0f92437..635af91 100644 --- a/lib/core/c++/chunk.hpp +++ b/lib/core/c++/chunk.hpp @@ -77,6 +77,11 @@ public: return values_.at(ind); } + /* + const data& neighbour_at(const data>& index) const { + } + */ + static constexpr auto get_dims(){ return data,Encode>{{Sides...}}; } diff --git a/lib/core/c++/hlbm.hpp b/lib/core/c++/hlbm.hpp index 7590cc2..ed49353 100644 --- a/lib/core/c++/hlbm.hpp +++ b/lib/core/c++/hlbm.hpp @@ -96,20 +96,55 @@ public: } }; -template +namespace impl { + +} + +template class component final { private: - typename saw::native_data_type::type relaxation_; - saw::data frequency_; + template + void apply_i(const saw::data& field, const saw::data& macros, const saw::data& part_groups, saw::data> index, saw::data time_step) const { + + } public: + /* + template + void apply_i() + */ template - void apply(const saw::data& field, const saw::data& macros, const saw::data& particles, saw::data> index, saw::data time_step) const { + void apply(const saw::data& field, const saw::data& macros, const saw::data& part_groups, saw::data> index, saw::data time_step) const { /// Figure out how to access the particle list // auto& p = particles.at(i); /// Iterate over the grid bounds // auto& grid = p.template get<"grid">(); + + auto& part_spheroid_group = part_groups.template get<0>(); + { + auto& parts = part_spheroid_group.template get<"particles">(); + auto parts_size = parts.size(); + + auto& pi = parts.at(index); + auto& pirb = pi.template get<"rigid_body">(); + auto& pirb_pos = pirb.template get<"position">(); + + saw::data> start; + saw::data> stop; + + /// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds + for(uint64_t i{0u}; i < Desc::D; ++i){ + + } + + iterator::apply([&](const auto& index){ + // ask for the d_k value here. + // For every value im iterating over I need sth + },start,stop); + + // Check + } } diff --git a/lib/core/c++/particle/aabb.hpp b/lib/core/c++/particle/aabb.hpp new file mode 100644 index 0000000..aec95ca --- /dev/null +++ b/lib/core/c++/particle/aabb.hpp @@ -0,0 +1,49 @@ +#pragma once + +#include "particle.hpp" + +namespace kel { +namespace lbm { +template +class particle_aabb final { +}; + +template::type radius> +class particle_aabb > > final { +public: + using Schema = sch::ParticleGroup>; + + using AABB = Struct< + Member"a">, + Member"b"> + >; +public: + static constexpr saw::data get(const saw::data& p_grp, const saw::data>& i, const saw::data>& meta){ + auto& parts = p_grp.template get<"particles">(); + + auto& pi = parts.at(i); + auto& pirb = pi.template get<"rigid_body">(); + auto& pirb_pos = pirb.template get<"position">(); + + saw::data aabb; + auto& a = aabb.template get<"a">(); + auto& b = aabb.template get<"b">(); + + saw::data> rad_d; + rad_d.at({}).set(radius); + + saw::data> lower; + saw::data> upper; + + for(uint64_t i{0u}; i < D; ++i){ + lower.at({{i}}) = pirb_pos.at({{i}}) >= rad_d.at({}) ? (pirb_pos.at({{i}}) - rad_d.at({})) : saw::data{0}; + a.at({i}) = lower.at({{i}}).template cast_to(); + upper.at({{i}}) = pirb_pos.at({{i}}) + rad_d.at({}); + b.at({i}) = (upper.at({{i}})+saw::data{1}).template cast_to() + } + + return aabb; + } +}; +} +} diff --git a/lib/core/c++/particle/porosity.hpp b/lib/core/c++/particle/porosity.hpp index aa1ce5b..39d9652 100644 --- a/lib/core/c++/particle/porosity.hpp +++ b/lib/core/c++/particle/porosity.hpp @@ -8,7 +8,9 @@ namespace lbm { template class particle_porosity { public: - saw::data> calculate(const saw::data<>& part_group, uint64_t p_i, const saw::data>& lbm_pos){ + + + saw::data> calculate(const saw::data>& part_group, uint64_t p_i, const saw::data>& lbm_pos){ auto& mask = part_group.template get<"mask">(); auto& particles = part_group.template get<"particles">(); @@ -17,7 +19,7 @@ public: auto& part_i_rb = part_i.template get<"rigid_body">(); auto& pirb = part_i_rb.template get<"position">(); - auto& dist = lbm_pos = lbm_pos - pirb; + auto dist = lbm_pos - pirb; // index 0 is at @@ -26,25 +28,58 @@ public: }; -template::type radius> -class particle_porosity> final { +template::type radius, typename saw::native_data_type::type eps> +class particle_porosity> final { public: - saw::data> calculate(const saw::data&, uint64_t i, const saw::data>& lbm_pos){ + saw::data> calculate(const saw::data > >& part_group, uint64_t i, const saw::data>& lbm_pos){ saw::data> por; por.at({}); - saw::data> dps_2; + auto& parts = part_group.template get<"particles">(); + auto& pi = parts.at({i}); + auto& pi_rb = pi.template get<"rigid_body">(); + + auto& pi_rb_pos = pi_rb.template get<"position">(); + + // Basically the queried position minus the center of the particle + auto dist = lbm_pos - pi_rb_pos; + + saw::data> dist_len; for(uint64_t i{0u}; i < D; ++i){ - auto& dps_i = lbm_pos.at({{i}}); - dps_2.at({}) = dps_i * dps_i; + auto& d_i = dist.at({{i}}); + dist_len.at({}) = d_i * d_i; } + dist_len.at({}).set(std::sqrt(dist_len.at({}).get()); - saw::data> rad_2; - rad_2.at({}).set(radius*radius); + saw::data> rad_d{radius}; + saw::data> eps_half{eps *0.5}; + + saw::data> rad_d_eps_p = rad_d + eps_half; + saw::data> rad_d_eps_n = rad_d - eps_half; + + // saw::data> rad_d_eps_p_2 = rad_d_eps_p * rad_d_eps_p; + // saw::data> rad_d_eps_n_2 = rad_d_eps_n * rad_d_eps_n; + // saw::data> rad_2; + // rad_2 = rad_d * rad_d; saw::data> inside; - if(dps_2.at({}).get() < rad_2.at({}).get()){ + if(dist_len.at({}) <= rad_d_eps_n.at({})){ inside.at({}).set(1); + return inside; + } + if(dist_len.at({}) > rad_d_eps_p.at({})){ + inside.at({}).set(0); + return inside + } + + /* + * cos^2 ( ||x-X(t)||_2 - (R-eps/2) ) + */ + { + typename saw::native_data_type::type inner = (std::pi / ( 2.0 * eps )) * (dist_len - rad_d_eps_n).at({}).get(); + auto cos_inner = std::cos(inner); + auto cos_i_2 = cos_inner * cos_inner; + inside.at({}).set(cos_i_2); } return inside; } -- cgit v1.2.3