diff options
Diffstat (limited to 'c++/particle')
| -rw-r--r-- | c++/particle/geometry/circle.hpp | 66 | ||||
| -rw-r--r-- | c++/particle/particle.hpp | 113 |
2 files changed, 0 insertions, 179 deletions
diff --git a/c++/particle/geometry/circle.hpp b/c++/particle/geometry/circle.hpp deleted file mode 100644 index 331f06d..0000000 --- a/c++/particle/geometry/circle.hpp +++ /dev/null @@ -1,66 +0,0 @@ -#pragma once - -#include "../particle.hpp" - -namespace kel { -namespace lbm { - -template<typename T> -class particle_circle_geometry { -private: -public: - particle_circle_geometry() - {} - - template<typename MT = T> - saw::data<sch::ParticleMask<MT,2>> generate_mask(uint64_t resolution, uint64_t boundary_nodes = 0) const { - saw::data<sch::ParticleMask<MT,2>> mask; - - auto& grid = mask.template get<"grid">(); - auto& com = mask.template get<"center_of_mass">(); - com = {}; - - //uint64_t rad_i = static_cast<uint64_t>(resolution * radius_.get())+1u; - uint64_t diameter_i = resolution; - uint64_t size = diameter_i + 2*boundary_nodes; - grid = {{{size,size}}}; - - saw::data<T> delta_x{static_cast<typename saw::native_data_type<T>::type>(2.0 / static_cast<double>(diameter_i))}; - saw::data<T> center = (saw::data<T>{1.0} + saw::data<T>{2.0} * saw::data<T>{static_cast<saw::native_data_type<T>::type>(boundary_nodes)/diameter_i}); - - for(uint64_t i = 0; i < size; ++i){ - for(uint64_t j = 0; j < size; ++j){ - grid.at({{i,j}}).set(0); - if(i >= boundary_nodes and j >= boundary_nodes and i < (diameter_i + boundary_nodes) and j < (diameter_i + boundary_nodes) ){ - saw::data<T> fi = saw::data<T>{static_cast<saw::native_data_type<T>::type>(0.5+i)} * delta_x - center; - saw::data<T> fj = saw::data<T>{static_cast<saw::native_data_type<T>::type>(0.5+j)} * delta_x - center; - - auto norm_f_ij = fi*fi + fj*fj; - if(norm_f_ij.get() <= 1){ - grid.at({{i,j}}).set(1); - } - } - } - } - - saw::data<T> total_mass{}; - iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){ - auto ind_vec = saw::math::vectorize_data(index).template cast_to<T>(); - for(uint64_t i{0u}; i < 2u; ++i){ - ind_vec.at({{i}}) = ind_vec.at({{i}}) * grid.at(index); - } - com = com + ind_vec; - - total_mass = total_mass + grid.at(index); - },{{0u,0u}}, grid.dims()); - - for(uint64_t i{0u}; i < 2u; ++i){ - com.at({{i}}) = com.at({{i}}) / total_mass; - } - - return mask; - } -}; - -} -} diff --git a/c++/particle/particle.hpp b/c++/particle/particle.hpp deleted file mode 100644 index 39aadfb..0000000 --- a/c++/particle/particle.hpp +++ /dev/null @@ -1,113 +0,0 @@ -#pragma once - -#include <forstio/codec/data.hpp> -#include <forstio/codec/data_math.hpp> -#include <forstio/codec/math.hpp> - -namespace kel { -namespace lbm { -namespace sch { -using namespace saw::schema; - -template<typename T, uint64_t D> -using ParticleRigidBody = Struct< - Member<Vector<T,D>, "position">, - Member<Vector<T,D>, "position_old">, - Member<T, "rotation">, - Member<T, "rotation_old">, - - Member<Vector<T,D>, "acceleration">, - Member<T, "rotational_acceleration"> ->; - -template<typename T, uint64_t D> -using ParticleMask = Struct< - Member<Array<T,D>, "grid">, - Member<Vector<T,D>, "center_of_mass"> ->; - -template<typename T, uint64_t D> -using Particle = Struct< - Member<ParticleRigidBody<T,D>, "rigid_body">, - Member<ParticleMask<T,D>, "mask">, - Member<T, "mass">, - Member<T, "size"> ->; -} - -template<typename T, uint64_t D, typename ParticleCollision = sch::ParticleMask<T,D> > -class particle_system { -private: - saw::data<sch::Array<sch::Particle<T,D>>> particles_; - - void verlet_step(saw::data<sch::Particle<T,D>>& particle, saw::data<T> time_step_delta){ - auto& body = particle.template get<"rigid_body">(); - - auto& pos = body.template get<"position">(); - auto& pos_old = body.template get<"position_old">(); - - // auto& rot = body.template get<"rotation">(); - auto& acc = body.template get<"acceleration">(); - - auto tsd_squared = time_step_delta * time_step_delta; - - saw::data<sch::Vector<T,D>> pos_new; - // Actual step - for(uint64_t i = 0u; i < D; ++i){ - pos_new.at({{i}}) = saw::data<T>{2.0} * pos.at({{i}}) - pos_old.at({{i}}) + acc.at({{i}}) * tsd_squared; - } - - pos_old = pos; - pos = pos_new; - } -public: - /** - * Add particle to this class and return an id representing this particle - */ - saw::error_or<saw::data<sch::UInt64>> add_particle(saw::data<sch::Particle<T,D>> particle__){ - auto size = particles_.size(); - auto eov = particles_.add(std::move(particle__)); - if(eov.is_error()){ - return std::move(eov.get_error()); - } - - return size; - } - - /* - saw::data<sch::Particle<T,D>>& get_particle(saw::data<sch::UInt64> id){ - } - */ - - void step(saw::data<T> time_step_delta){ - for(saw::data<sch::UInt64> i{0u}; i < particles_.size(); ++i){ - verlet_step(particles_.at(i), time_step_delta); - } - } - - template<typename LbmLattice> - void update_particle_border(saw::data<LbmLattice>& latt){ - (void) latt; - for(auto& iter : particles_){ - auto& par = iter; - - auto& body = par.template get<"rigid_body">(); - auto& size = par.template get<"size">(); - - - } - } - - saw::data<sch::UInt64> size() const { - return particles_.size(); - } - - /** - * Mostly meant for unforeseen use cases. - */ - saw::data<sch::Particle<T,D>>& at(saw::data<sch::UInt64> i){ - return particles_.at(i); - } -}; -} -} |
