1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
#pragma once
#include <forstio/codec/data.hpp>
#include <forstio/codec/data_math.hpp>
#include <forstio/codec/math.hpp>
namespace kel {
namespace lbm {
namespace coll {
struct Spheroid{};
}
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>
using ParticleCollisionSpheroid = Struct<
Member<T, "radius">
>;
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>
constexpr auto verlet_step_lambda = [](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;
};
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);
}
};
}
}
|
