-
Notifications
You must be signed in to change notification settings - Fork 0
/
singularstate.cpp
143 lines (120 loc) · 4.42 KB
/
singularstate.cpp
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
#include "singularstate.h"
#include <QtCore/qmath.h>
#include <random>
namespace Simulation {
SingularState::SingularState() :
particles()
{
}
void SingularState::add_particle(QPointF point, const VelocityMap *velocity_map)
{
Particle new_p;
new_p.position = point;
new_p.velocity = velocity_map->Get_Velocity_At(point.x(), point.y());
particles.append(new_p);
}
void SingularState::remove_particles(QList<int> positions)
{
std::sort(positions.begin(), positions.end());
positions.append(particles.length());
int ifrom = 0, ito = 0;
for (int i: positions)
{
while (ifrom < i)
{
std::swap(particles[ifrom++], particles[ito++]);
}
ifrom++;
}
particles.remove(ito, (ifrom-1)-ito);
}
QList<int> SingularState::find_enclosed(QPolygonF shape)
{
QList<int> res;
int l = particles.size();
for (int i = 0; i < l; i++) {
if (shape.containsPoint(particles[i].position, Qt::OddEvenFill))
res.append(i);
}
return res;
}
void SingularState::generate_enclosed(QPolygonF shape, qreal density, const VelocityMap *velocity_map)
{
//Calculating the amount of points needed.
qreal area = 0;
QPointF s = shape[0];
for (int i = 2; i < shape.size(); i++) {
QPointF v1 = shape[i-1] - s, v2 = shape[i] - s;
area += v1.x() * v2.y() - v1.y() * v2.x();
}
area = qFabs(area);
int number_of_points = qRound(area * density);
QRectF rect = shape.boundingRect();
qreal x1, y1, x2, y2;
rect.getCoords(&x1, &y1, &x2, &y2);
// The following is a bit cheesy, but indubitably simple and fast enough for real-life purposes.
std::random_device rd;
std::default_random_engine engine(rd());
std::uniform_real_distribution<qreal> randx(x1, std::nextafter(x2, std::numeric_limits<qreal>::max()));
std::uniform_real_distribution<qreal> randy(y1, std::nextafter(y2, std::numeric_limits<qreal>::max()));
while (number_of_points)
{
QPointF p(randx(engine), randy(engine));
if (shape.containsPoint(p, Qt::OddEvenFill))
{
number_of_points--;
add_particle(p, velocity_map);
}
}
}
void SingularState::update_from(const SingularState *ancestor, const VelocityMap *velocity_map)
{
int old_size = particles.size(), new_size = ancestor->particles.size();
// here, assuming new_size >= old_size
particles.resize(new_size);
std::copy(ancestor->particles.begin() + old_size, ancestor->particles.end(), particles.begin() + old_size);
recalculate(old_size, new_size, velocity_map);
}
void SingularState::full_update(const SingularState *ancestor, const VelocityMap *velocity_map)
{
particles.resize(ancestor->particles.size());
std::copy(ancestor->particles.begin(), ancestor->particles.end(), particles.begin());
recalculate(0, particles.size(), velocity_map);
}
SingularState SingularState::successor(const VelocityMap *velocity_map)
{
SingularState result;
result.full_update(this, velocity_map);
return result;
}
void SingularState::recalculate(int from, int to, const VelocityMap *velocity_map)
{
for (int i = from; i < to; i++) {
qreal h = STEP;
Particle &p = particles[i];
QPointF pos[4], vel[4], force[4];
pos[0] = p.position; vel[0] = p.velocity; force[0] = velocity_map->Calculate_Force(p.position.x(), p.position.y(), p.velocity);
pos[1] = pos[0] + h/2 * vel[0]; vel[1] = vel[0] + h/2 * force[0]; force[1] = velocity_map->Calculate_Force(pos[1].x(), pos[1].y(), vel[1]);
pos[2] = pos[0] + h/2 * vel[1]; vel[1] = vel[0] + h/2 * force[1]; force[2] = velocity_map->Calculate_Force(pos[2].x(), pos[2].y(), vel[2]);
pos[3] = pos[0] + h * vel[2]; vel[1] = vel[0] + h * force[2]; force[3] = velocity_map->Calculate_Force(pos[3].x(), pos[3].y(), vel[3]);
p.position += h * ((vel[0] + vel[3])/ 6 + (vel[1] + vel[2]) / 3);
p.velocity += h * ((force[0] + force[3]) / 6 + (force[1] + force[2]) / 3);
}
}
QDataStream &operator>>(QDataStream &input, Particle &self)
{
return input >> self.position >> self.velocity;
}
QDataStream &operator<<(QDataStream &output, const Particle &self)
{
return output << self.position << self.velocity;
}
QDataStream& operator>>(QDataStream &input, SingularState &self)
{
return input >> self.particles;
}
QDataStream& operator<<(QDataStream &output, const SingularState &self)
{
return output << self.particles;
}
}