GCC Code Coverage Report

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