GCC Code Coverage Report

Directory: ./
File: src/system.cpp
Date: 2025-11-05 16:49:12
Exec Total Coverage
Lines: 185 230 80.4%
Functions: 16 18 88.9%
Branches: 171 299 57.2%
Line Branch Exec Source
1 #include <fstream>
2 #include <iostream>
3 #include <limits>
4 #include <sstream>
5
6 #include "analyzers.hpp"
7 #include "interaction.hpp"
8 #include "multiindex.hpp"
9 #include "particle.hpp"
10 #include "rng.hpp"
11 #include "system.hpp"
12
13 75 System::System(std::vector<Particle> &particles, const YAML::Node &config, double t,
14 75 const YAML::Node &addParticlesConfig)
15
1/2
✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
 150 : particles(particles), t(t), L(YAML::parse<Vec>(config, "L")),
16
1/2
✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
 150 T(YAML::parse<double>(config, "T", std::numeric_limits<double>::quiet_NaN())),
17
1/2
✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
 150 mu(YAML::parse<double>(config, "mu", std::numeric_limits<double>::quiet_NaN())), E(0.),
18
2/4
✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 75 times.
✗ Branch 7 not taken.
 75 interaction(InteractionFactory::create(this, config["interaction"])) {
19
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 std::clog << "\nSystem:";
20
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 std::clog << "\n\tL:";
21
2/2
✓ Branch 0 taken 225 times.
✓ Branch 1 taken 75 times.
 300 for (int i = 0; i < DIM; ++i) {
22
2/4
✓ Branch 1 taken 225 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 225 times.
✗ Branch 6 not taken.
 225 std::clog << " " << this->L[i];
23 }
24
1/2
✓ Branch 0 taken 75 times.
✗ Branch 1 not taken.
 75 if (std::isfinite(T)) {
25
2/4
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
 75 std::clog << "\n\tT: " << this->T;
26 }
27
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 25 times.
 75 if (std::isfinite(mu)) {
28
2/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
 50 std::clog << "\n\tmu: " << this->mu;
29 }
30
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 std::clog << "\nParticles:";
31
3/6
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 75 times.
✗ Branch 8 not taken.
 75 std::clog << "\n\t" << particles.size() << " already in the system";
32
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 initParticles(addParticlesConfig);
33
5/10
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 75 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 75 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 75 times.
✗ Branch 13 not taken.
 150 std::clog << "\n\ttotal: " << particles.size() << " (N/V = " << particles.size() / L.prod()
34
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 << ")";
35
2/4
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 75 times.
 75 if (!allParticlesInside()) {
36 ✗ throw std::runtime_error(
37 ✗ "Could not initialize the system because some particles would be outside of the box.");
38 }
39
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 25 times.
 75 if (particles.size() == 0) {
40 50 std::clog << "\nWARNING: Initialized system with zero particles. Ensure that this is intended, "
41
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 "otherwise check the particle initialization in the config!";
42 }
43 75 }
44
45
1/2
✓ Branch 0 taken 73 times.
✗ Branch 1 not taken.
 73 System::~System() { delete interaction; }
46
47 ✗ void System::ensureAnalyzer(const std::string &key) {
48 ✗ if (!analyzers[key]) {
49 ✗ analyzers[key] = AnalyzerFactory::create(this, key);
50 }
51 }
52
53 ✗ void System::updateAnalyzers() {
54 ✗ for (auto &[_, analyzer] : analyzers) {
55 ✗ analyzer->update();
56 }
57 }
58
59 10 int System::getNumSpecies() {
60 10 int maxSpecies = 0;
61
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 10 times.
 12 for (auto &p : particles) {
62
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 maxSpecies = std::max(maxSpecies, p.species);
63 }
64 10 return maxSpecies + 1;
65 }
66
67 4885 bool System::removeLastParticleIfHighEnergy(double threshold) {
68 4885 bool wasRemoved = false;
69 4885 try {
70 4885 double Eparticle = 0.;
71
1/2
✓ Branch 1 taken 4885 times.
✗ Branch 2 not taken.
 4885 Eparticle = updateParticleState(particles.back(), ParticleUpdateType::INSERT);
72
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4885 times.
 4885 if (Eparticle >= threshold) {
73 ✗ std::clog << "\n\tParticle insertion skipped due to energy threshold (E=" << Eparticle
74 ✗ << " >= threshold=" << threshold << ")";
75 ✗ updateParticleMeta(particles.back(), ParticleUpdateType::DELETE);
76 ✗ particles.pop_back();
77 ✗ wasRemoved = true;
78 }
79 ✗ } catch (const std::runtime_error &e) {
80 ✗ std::clog << "\nWARNING: Could not check energy threshold of particle insertion because "
81 ✗ "interaction does not support single particle updates";
82 }
83 4885 return wasRemoved;
84 }
85
86 78 void System::initParticles(const YAML::Node &addParticlesConfig) {
87
1/2
✓ Branch 0 taken 78 times.
✗ Branch 1 not taken.
 78 if (addParticlesConfig) {
88
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 77 times.
 78 if (addParticlesConfig.IsSequence()) {
89
4/6
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✓ Branch 9 taken 1 times.
 3 for (auto &c : addParticlesConfig) {
90
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 initParticlesItem(c);
91
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 3 }
92 } else {
93 77 initParticlesItem(addParticlesConfig);
94 }
95 }
96 78 }
97
98 78 void System::initParticlesItem(const YAML::Node &addParticlesConfig) {
99 78 int nBefore = particles.size();
100 78 double defaultEThreshold = 1e4;
101
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 78 times.
 78 if (addParticlesConfig["file"]) {
102 ✗ auto &config = addParticlesConfig["file"];
103 ✗ std::string filename = YAML::parse<std::string>(config, "filename");
104 ✗ std::ifstream particlesInitStream(filename, std::ios_base::in);
105 ✗ if (particlesInitStream.fail()) {
106 ✗ throw std::invalid_argument("Could not open initialization file " + filename);
107 }
108 ✗ double EThreshold = YAML::parse<double>(config, "EThreshold", defaultEThreshold);
109 ✗ initParticlesFromFile(particlesInitStream, EThreshold);
110 ✗ std::clog << "\n\tAdded " << particles.size() - nBefore << " from file " << filename;
111 }
112
2/2
✓ Branch 2 taken 28 times.
✓ Branch 3 taken 50 times.
 78 if (addParticlesConfig["lattice"]) {
113 28 auto &config = addParticlesConfig["lattice"];
114
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
 56 std::string cell = YAML::parse<std::string>(config, "cell", "sc");
115
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
 56 Vec spacing = YAML::parse<Vec>(config, "spacing", 1.);
116
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
 56 Vec margin = YAML::parse<Vec>(config, "margin", 0.5);
117
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
 56 Vec offset = YAML::parse<Vec>(config, "offset", 0.);
118
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
 56 int NTarget = YAML::parse<int>(config, "N", 0);
119
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
 56 double EThreshold = YAML::parse<double>(config, "EThreshold", defaultEThreshold);
120
1/2
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
 28 YAML::Node particleProperties;
121
2/4
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
 28 if (config["particleProperties"]) {
122
2/4
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
 28 particleProperties = config["particleProperties"];
123 }
124
1/2
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
 28 initParticlesOnLattice(cell, spacing, margin, offset, particleProperties, NTarget, EThreshold);
125
5/10
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 28 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 28 times.
✗ Branch 14 not taken.
 28 std::clog << "\n\tAdded " << particles.size() - nBefore << " on " << cell << " lattice";
126 56 }
127 78 }
128
129 16 void System::initParticlesFromFile(std::istream &stream, double EThreshold) {
130 // Not very robust, I know
131 16 std::string line;
132
3/4
✓ Branch 1 taken 171 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 155 times.
✓ Branch 4 taken 16 times.
 171 while (std::getline(stream, line)) {
133
4/4
✓ Branch 0 taken 153 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 152 times.
 155 if (line.empty() || line.front() == '#') {
134 3 continue;
135 }
136
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 152 std::istringstream iss(line);
137
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 152 iss.exceptions(std::ios::failbit);
138 152 Vec r;
139 152 Vec v = Vec::Constant(0.);
140 152 Vec u = Vec::Constant(0.);
141 152 Real sigma = 1;
142 152 Real m = 1;
143 152 Real gamma = 1;
144 152 int species = 0;
145
2/2
✓ Branch 0 taken 456 times.
✓ Branch 1 taken 152 times.
 608 for (int i = 0; i < DIM; ++i) {
146 456 Real ri;
147
1/2
✓ Branch 1 taken 456 times.
✗ Branch 2 not taken.
 456 iss >> ri;
148
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 456 times.
 456 if (ri >= L[i]) {
149 ✗ throw std::invalid_argument("Particle " + std::to_string(i) +
150 ✗ " is not inside simulation box. Aborting");
151 }
152 456 r[i] = ri;
153 }
154 try {
155
2/2
✓ Branch 0 taken 456 times.
✓ Branch 1 taken 152 times.
 608 for (int i = 0; i < DIM; ++i) {
156
1/2
✓ Branch 2 taken 456 times.
✗ Branch 3 not taken.
 456 iss >> v[i];
157 }
158
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 152 iss >> sigma;
159
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 152 iss >> m;
160
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 152 iss >> gamma;
161
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 152 times.
 152 iss >> species;
162
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 152 times.
 152 } catch (const std::istream::failure &e) {
163 // ignore optional properties
164 152 }
165
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 152 particles.emplace_back(r, v, u, sigma, m, gamma, species);
166
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 152 removeLastParticleIfHighEnergy(EThreshold);
167 152 }
168 16 }
169
170 28 static void removeAverageVelocity(std::vector<Particle> &particles) {
171 28 Vec vAvg = Vec::Constant(0.);
172
2/2
✓ Branch 0 taken 4733 times.
✓ Branch 1 taken 28 times.
 4761 for (const Particle &p : particles) {
173 4733 vAvg += p.v;
174 }
175 28 vAvg /= particles.size();
176
2/2
✓ Branch 0 taken 4733 times.
✓ Branch 1 taken 28 times.
 4761 for (Particle &p : particles) {
177 4733 p.v -= vAvg;
178 }
179 28 }
180
181 28 void System::initParticlesOnLattice(const std::string &cell, Vec spacing, const Vec &margin,
182 const Vec &offset, const YAML::Node &particleProperties,
183 int NTarget, double EThreshold) {
184
2/4
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 28 times.
 28 if ((spacing.array() <= 0.).any() || (spacing.array() == 1.).all()) {
185 ✗ spacing = Vec::Constant(1. + 1e-10); // Safety margin due to floating point errors
186 }
187
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 27 times.
 28 Vec a = spacing;
188
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 27 times.
 28 if (cell == "bcc") {
189 1 a *= 2. / std::sqrt(3.);
190
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 26 times.
 27 } else if (cell == "fcc") {
191 1 a *= std::sqrt(2.);
192 }
193
2/2
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 2 times.
 28 std::vector<Vec> bravais;
194 28 if constexpr (DIM == 2) {
195 if (cell == "sc") {
196 bravais.push_back(Vec::Constant(0.));
197 } else {
198 throw std::invalid_argument("Grid type " + cell + " not implemented!");
199 }
200 }
201 28 if constexpr (DIM == 3) {
202
2/2
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 2 times.
 28 if (cell == "sc") {
203
1/2
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
 26 bravais.push_back(Vec::Constant(0.));
204
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 } else if (cell == "bcc") {
205
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 bravais.push_back(Vec::Constant(0.));
206
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 bravais.push_back(Vec{0.5, 0.5, 0.5});
207
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 } else if (cell == "fcc") {
208
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 bravais.push_back(Vec::Constant(0.));
209
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 bravais.push_back(Vec{0., 0.5, 0.5});
210
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 bravais.push_back(Vec{0.5, 0., 0.5});
211
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 bravais.push_back(Vec{0.5, 0.5, 0.});
212 } else {
213 ✗ throw std::invalid_argument("Grid type " + cell + " not implemented!");
214 }
215 }
216
217
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
 56 Real sigma = YAML::parse<Real>(particleProperties, "sigma", 1.);
218
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
 56 Real m = YAML::parse<Real>(particleProperties, "m", 1.);
219
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 28 times.
✗ Branch 9 not taken.
 56 Real gamma = YAML::parse<Real>(particleProperties, "gamma", 1.);
220
2/4
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
 56 int species = YAML::parse<int>(particleProperties, "species", 0);
221 28 bool vMaxwell = false;
222 28 bool uRandom = false;
223 28 Vec v, u = Vec::Zero();
224 28 try {
225
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 28 times.
 84 if (YAML::parse<std::string>(particleProperties, "v", "") == "maxwell") {
226 ✗ vMaxwell = true;
227 }
228 ✗ } catch (const YAML::BadConversion &e) {
229 ✗ v = YAML::parse<Vec>(particleProperties, "v", 0.);
230 }
231 28 try {
232
3/6
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 28 times.
 84 if (YAML::parse<std::string>(particleProperties, "u", "") == "random") {
233 ✗ uRandom = true;
234 }
235 ✗ } catch (const YAML::BadConversion &e) {
236 ✗ u = YAML::parse<Vec>(particleProperties, "u", 0.);
237 }
238
239
1/2
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
 28 RNG rng;
240 28 std::gamma_distribution<double> maxwellDist(1.5, T);
241
242 28 size_t NSpecies = 0;
243 28 ArrayI bounds = (L.array() / a.array() + 1).cast<int>();
244
2/2
✓ Branch 1 taken 10174 times.
✓ Branch 2 taken 28 times.
 10202 for (ArrayI multi : Multiindex<DIM>(bounds)) {
245 10174 Vec rc = multi.cast<Real>() * a.array();
246
2/2
✓ Branch 0 taken 10342 times.
✓ Branch 1 taken 10174 times.
 20516 for (auto &b : bravais) {
247 10342 Vec r = rc + (b.array() * a.array()).matrix() + offset;
248
2/2
✓ Branch 0 taken 4733 times.
✓ Branch 1 taken 5609 times.
 10342 if ((r.array() >= margin.array() && r.array() < L.array() - margin.array()).all()) {
249
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4733 times.
 4733 if (vMaxwell) {
250 ✗ v = rng.vecNormal().normalized();
251 ✗ v *= maxwellDist(rng.engine);
252 }
253
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4733 times.
 4733 if (uRandom) {
254 ✗ do {
255 ✗ u = rng.vecNormal();
256 ✗ } while (u.squaredNorm() - 1. > 1e-5);
257 }
258
1/2
✓ Branch 1 taken 4733 times.
✗ Branch 2 not taken.
 4733 u.normalize();
259
1/2
✓ Branch 1 taken 4733 times.
✗ Branch 2 not taken.
 4733 particles.emplace_back(r, v, u, sigma, m, gamma, species);
260
2/4
✓ Branch 1 taken 4733 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4733 times.
✗ Branch 4 not taken.
 4733 if (!removeLastParticleIfHighEnergy(EThreshold)) {
261 4733 ++NSpecies;
262 }
263
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 4733 times.
✗ Branch 3 not taken.
 4733 if (NTarget && NSpecies == (size_t)NTarget) {
264 ✗ removeAverageVelocity(particles);
265 ✗ return;
266 }
267 }
268 }
269 }
270 28 removeAverageVelocity(particles);
271
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 28 times.
✗ Branch 3 not taken.
 28 if (NTarget && NSpecies < (size_t)NTarget) {
272 ✗ throw std::runtime_error("Could not fit " + std::to_string(NTarget) + " particles of species" +
273 ✗ std::to_string(species) + " into the box (only got to " +
274 ✗ std::to_string(NSpecies) + ").");
275 }
276 28 }
277
278 10 void System::updateParticleMeta(Particle &p, ParticleUpdateType particleUpdateType) {
279 10 interaction->updateParticleMeta(p, particleUpdateType);
280 10 }
281
282 5102 double System::updateParticleState(Particle &p, ParticleUpdateType particleUpdateType) {
283 5102 p.Fext.setZero();
284 5102 p.Eext = 0.;
285 5102 p.Fint.setZero();
286 5102 return interaction->updateParticleState(p, particleUpdateType);
287 }
288
289 4789 void System::updateState() {
290 4789 E = 0.;
291
2/2
✓ Branch 0 taken 18928 times.
✓ Branch 1 taken 4789 times.
 23717 for (auto &p : particles) {
292 18928 p.Fext.setZero();
293 18928 p.Eext = 0.;
294 18928 p.Fint.setZero();
295 }
296 4789 interaction->updateState();
297 4789 }
298
299 14944 bool System::particleInside(const Particle &p) {
300
2/2
✓ Branch 1 taken 3608 times.
✓ Branch 2 taken 11336 times.
 14944 if ((p.r.array() >= L.array() || p.r.array() < 0).any()) {
301 3608 return false;
302 }
303 return true;
304 }
305
306 78 bool System::allParticlesInside() {
307
2/2
✓ Branch 0 taken 4733 times.
✓ Branch 1 taken 78 times.
 4811 for (auto &p : particles) {
308
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4733 times.
 4733 if (!particleInside(p)) {
309 ✗ return false;
310 }
311 }
312 78 return true;
313 }
314
315 6603 void System::applyPeriodic(Particle &p) {
316
2/2
✓ Branch 1 taken 3608 times.
✓ Branch 2 taken 6603 times.
 10211 while (!particleInside(p)) {
317 3608 Vec shift = (p.r.array() / L.array()).floor();
318 3608 p.box += shift.cast<int>().array();
319 7216 p.r.array() -= shift.array() * L.array();
320 }
321 6603 }
322
323 24 void System::applyPeriodic() {
324 24 #pragma omp parallel for
325 for (auto &p : particles) {
326 applyPeriodic(p);
327 }
328 24 }
329