diff options
Diffstat (limited to 'lib/core/c++')
| -rw-r--r-- | lib/core/c++/boundary.hpp | 8 | ||||
| -rw-r--r-- | lib/core/c++/particle/particle.hpp | 101 |
2 files changed, 85 insertions, 24 deletions
diff --git a/lib/core/c++/boundary.hpp b/lib/core/c++/boundary.hpp index 4dbbdf8..01ae7b5 100644 --- a/lib/core/c++/boundary.hpp +++ b/lib/core/c++/boundary.hpp @@ -211,12 +211,12 @@ public: if constexpr (Dir) { dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data<FP>{2.0 / 3.0} * rho_vel_x; - dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); - dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); + dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); + dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); }else if constexpr (not Dir){ dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data<FP>{2.0 / 3.0} * rho_vel_x; - dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); - dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); + dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); + dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); } } }; diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp index f8104fd..c0d115f 100644 --- a/lib/core/c++/particle/particle.hpp +++ b/lib/core/c++/particle/particle.hpp @@ -1,8 +1,10 @@ #pragma once -#include <forstio/codec/data.hpp> -#include <forstio/codec/data_math.hpp> #include <forstio/codec/math.hpp> +#include <forstio/codec/data_math.hpp> +#include <forstio/codec/data.hpp> + +#include "../iterator.hpp" namespace kel { namespace lbm { @@ -43,19 +45,64 @@ using ParticleCollisionSpheroid = Struct< Member<Scalar<T>, "radius"> >; -tem - -template<typename T, uint64_t D, typename CollisionType = ParticleCollisionSpheroid<T>> +template<typename T, uint64_t D> using Particle = Struct< - Member<ParticleRigidBody<T,D>, "rigid_body">, - Member<CollisionType, "collision">, + Member<ParticleRigidBody<T,D>, "rigid_body"> // Problem is that dynamic data would two layered // Member<Array<Float64,D>, "mask">, - Member<UInt64, "mask_id">, - Member<Scalar<T>, "mass"> >; + +template<typename T, uint64_t D, typename CollisionType = ParticleCollisionSpheroid<T>> +using ParticleGroup = Struct< + Member<Array<T,D>, "mask">, + Member<CollisionType, "collision">, + Member<Scalar<T>, "mass">, + Member<Array<Particle<T,D>>, "particles"> +>; +} + +template<typename T, uint64_t D> +saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T>>> create_spheroid_particle_group( + saw::data<sch::Scalar<T>> rad_p, + saw::data<sch::Scalar<T>> density_p, + const saw::data<sch::UInt64>& mask_resolution +){ + saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T>>> part; + + auto& mask = part.template get<"mask">(); + auto& collision = part.template get<"collision">(); + auto& mass = part.template get<"mass">(); + + if constexpr ( D == 1u ){ + saw::data<sch::Scalar<T>> c; + c.at({}).set(2.0); + mass = rad_p * c * density_p; + } else if constexpr ( D == 2u){ + saw::data<sch::Scalar<T>> pi; + pi.at({}).set(3.14159); + mass = rad_p * rad_p * pi * density_p; + } else if constexpr ( D == 3u ){ + saw::data<sch::Scalar<T>> c; + c.at({}).set(3.14159 * 4.0 / 3.0); + mass = rad_p * rad_p * rad_p * c * density_p; + } else { + static_assert(D == 0u or D > 3u, "Dimensions only supported for Dim 1,2 & 3."); + } + + saw::data<sch::FixedArray<sch::UInt64,D>> mask_dims; + for(uint64_t i = 0u; i < D; ++i){ + mask_dims.at({i}) = mask_resolution; + } + + mask = {mask_dims}; + iterator<D>::apply([&](const auto& index){ + + },{},mask_dims); + + return part; } +/* template<typename T, uint64_t D> saw::data<sch::Particle<T,D, sch::ParticleCollisionSpheroid<T>>> create_spheroid_particle( saw::data<sch::Vector<T,D>> pos_p, @@ -106,6 +153,7 @@ saw::data<sch::Particle<T,D, sch::ParticleCollisionSpheroid<T>>> create_spheroid return part; } +*/ template<typename T,uint64_t D> constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle, saw::data<sch::Scalar<T>> time_step_delta){ @@ -142,18 +190,20 @@ constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle, }; template<typename T, uint64_t D> -constexpr auto handle_collision = [](saw::data<sch::Particle<T,D>>& left, saw::data<sch::Particle<T,D>>& right, saw::data<sch::Vector<T,D>> unit_pos_rel, saw::data<sch::Vector<T,D>> vel_rel, saw::data<sch::Scalar<T>> d_t) -{ +constexpr auto handle_collision = []( + saw::data<sch::Particle<T,D>>& left, const saw::data<sch::Scalar<T>>& mass_l, + saw::data<sch::Particle<T,D>>& right, const saw::data<sch::Scalar<T>>& mass_r, + saw::data<sch::Vector<T,D>> unit_pos_rel, saw::data<sch::Vector<T,D>> vel_rel, + saw::data<sch::Scalar<T>> d_t +){ auto& rb_l = left.template get<"rigid_body">(); auto& pos_l = rb_l.template get<"position">(); auto& pos_old_l = rb_l.template get<"position_old">(); - auto& mass_l = left.template get<"mass">(); auto vel_l = (pos_l-pos_old_l)/d_t; auto& rb_r = right.template get<"rigid_body">(); auto& pos_r = rb_r.template get<"position">(); auto& pos_old_r = rb_r.template get<"position_old">(); - auto& mass_r = left.template get<"mass">(); auto vel_r = (pos_r-pos_old_r)/d_t; auto vel_pos_rel_dot = saw::math::dot(unit_pos_rel,vel_rel); @@ -165,10 +215,16 @@ constexpr auto handle_collision = [](saw::data<sch::Particle<T,D>>& left, saw::d }; -template<typename T, uint64_t D> -constexpr auto broadphase_collision_distance_squared = [](saw::data<sch::Particle<T,D>>& left, saw::data<sch::Particle<T,D>>& right) -> std::pair<bool,saw::data<sch::Scalar<T>>>{ - auto rad_l = left.template get<"collision">().template get<"radius">(); - auto rad_r = right.template get<"collision">().template get<"radius">(); +template<typename T, uint64_t D, typename Collision> +constexpr auto broadphase_collision_distance_squared = []( + saw::data<sch::Particle<T,D>>& left, + const saw::data<Collision>& coll_l, + saw::data<sch::Particle<T,D>>& right, + const saw::data<Collision>& coll_r +) -> std::pair<bool,saw::data<sch::Scalar<T>>>{ + + auto rad_l = coll_l.template get<"radius">(); + auto rad_r = coll_r.template get<"radius">(); auto& rb_l = left.template get<"rigid_body">(); auto& rb_r = right.template get<"rigid_body">(); @@ -190,9 +246,14 @@ constexpr auto broadphase_collision_distance_squared = [](saw::data<sch::Particl * * */ -template<typename T, uint64_t D> -constexpr auto broadphase_collision_check = [](saw::data<sch::Particle<T,D>>& left, saw::data<sch::Particle<T,D>>& right) -> bool{ - return broadphase_collision_distance_squared<T,D>(left,right).first; +template<typename T, uint64_t D, typename Collision> +constexpr auto broadphase_collision_check = []( + saw::data<sch::Particle<T,D>>& left, + const saw::data<Collision>& coll_l, + saw::data<sch::Particle<T,D>>& right, + const saw::data<Collision>& coll_r +) -> bool{ + return broadphase_collision_distance_squared<T,D,Collision>(left,coll_l,right,coll_r).first; }; } } |