GCC Code Coverage Report

Directory: ./
File: src/cells.cpp
Date: 2025-07-07 11:13:23
Exec Total Coverage
Lines: 124 128 96.9%
Functions: 5 5 100.0%
Branches: 123 168 73.2%
Line Branch Exec Source
1 #include <iostream>
2 #include <stdexcept>
3
4 #include "cells.hpp"
5
6 88 static int getNumColors(const std::string &mode) {
7 88 int colors;
8
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 88 times.
 88 if (mode == "c08") {
9 colors = 8;
10 ✗ } else if (mode == "c04") {
11 colors = 4;
12 } else {
13 ✗ throw std::invalid_argument("Cell color scheme " + mode + " not implemented");
14 }
15 88 return colors;
16 }
17
18 88 Cells::Cells(const Vec &cellL, const ArrayI &count, const std::string &colorScheme, bool triples)
19
2/4
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 88 times.
✗ Branch 7 not taken.
 88 : cellLinv(1. / cellL.array()), multiindex(count), colorScheme(colorScheme) {
20
1/2
✓ Branch 0 taken 88 times.
✗ Branch 1 not taken.
 88 if (this->colorScheme.empty()) {
21 // Set default color scheme if none given
22 88 if constexpr (DIM == 3) {
23
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 this->colorScheme = "c08";
24 }
25 if constexpr (DIM == 2) {
26 this->colorScheme = "c04";
27 }
28 }
29
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 this->cells.resize(multiindex.linsize);
30
2/2
✓ Branch 0 taken 7550 times.
✓ Branch 1 taken 88 times.
 7638 for (int i = 0; i < multiindex.linsize; ++i) {
31
1/2
✓ Branch 1 taken 7550 times.
✗ Branch 2 not taken.
 7550 this->cells[i].particleIndices = {};
32
3/6
✓ Branch 1 taken 7550 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7550 times.
✗ Branch 5 not taken.
✓ Branch 10 taken 7550 times.
✗ Branch 11 not taken.
 15100 this->cells[i].pairsAll = getPairs(i, "all");
33
1/2
✓ Branch 1 taken 7550 times.
✗ Branch 2 not taken.
 7550 this->cells[i].pairsNewton3 = getPairs(i, this->colorScheme);
34
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 7496 times.
 7550 if (triples) {
35
3/6
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
✓ Branch 10 taken 54 times.
✗ Branch 11 not taken.
 108 this->cells[i].triplesAll = getTriples(i, "all");
36
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 this->cells[i].triplesNewton3 = getTriples(i, this->colorScheme);
37 }
38 }
39
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 int colors = getNumColors(this->colorScheme);
40
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 this->sameColorIndices.resize(colors);
41
2/2
✓ Branch 0 taken 704 times.
✓ Branch 1 taken 88 times.
 792 for (int c = 0; c < colors; c++) {
42
1/2
✓ Branch 1 taken 704 times.
✗ Branch 2 not taken.
 1408 this->sameColorIndices.push_back(getSameColorIndices(c, this->colorScheme));
43 }
44
2/4
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
 88 std::clog << "\n\tcell color scheme: " << this->colorScheme;
45
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 std::clog << "\n\tcell count:";
46
2/2
✓ Branch 0 taken 264 times.
✓ Branch 1 taken 88 times.
 352 for (int d = 0; d < DIM; ++d) {
47
2/4
✓ Branch 1 taken 264 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 264 times.
✗ Branch 6 not taken.
 264 std::clog << " " << this->multiindex.shape[d];
48 }
49
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 std::clog << "\n\tcell size:";
50
2/2
✓ Branch 0 taken 264 times.
✓ Branch 1 taken 88 times.
 352 for (int d = 0; d < DIM; ++d) {
51
2/4
✓ Branch 1 taken 264 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 264 times.
✗ Branch 6 not taken.
 264 std::clog << " " << cellL[d];
52 }
53 88 }
54
55 15103 std::vector<std::pair<int, int>> Cells::getPairs(int cellIndex, const std::string &mode) {
56 15103 ArrayI cellMulti = multiindex.lin2multi(cellIndex);
57
2/2
✓ Branch 0 taken 7550 times.
✓ Branch 1 taken 7553 times.
 15103 std::vector<std::pair<ArrayI, ArrayI>> shiftPairs;
58
2/2
✓ Branch 0 taken 7550 times.
✓ Branch 1 taken 7553 times.
 15103 if (mode == "all") {
59 7550 ArrayI shift1 = ArrayI::Zero();
60
2/2
✓ Branch 2 taken 203850 times.
✓ Branch 3 taken 7550 times.
 211400 for (ArrayI shift2 : Multiindex<DIM>(ArrayI(3))) {
61 203850 shift2 -= ArrayI(1);
62
1/2
✓ Branch 1 taken 203850 times.
✗ Branch 2 not taken.
 203850 shiftPairs.push_back({shift1, shift2});
63 }
64
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 7553 times.
✗ Branch 3 not taken.
 7553 } else if (mode == "c08" || mode == "c04") {
65 7553 Multiindex<DIM> range(ArrayI(2));
66
2/2
✓ Branch 0 taken 60424 times.
✓ Branch 1 taken 7553 times.
 67977 for (ArrayI shift1 : range) {
67
2/2
✓ Branch 0 taken 483392 times.
✓ Branch 1 taken 60424 times.
 543816 for (ArrayI shift2 : range) {
68
2/2
✓ Branch 1 taken 279461 times.
✓ Branch 2 taken 203931 times.
 483392 if ((shift1 * shift2 == 1).any()) {
69 279461 continue;
70 }
71
1/2
✓ Branch 1 taken 203931 times.
✗ Branch 2 not taken.
 203931 shiftPairs.push_back({shift1, shift2});
72 }
73 }
74 } else if (mode == "c04") {
75 throw std::invalid_argument("Cell mode " + mode + " not implemented");
76 }
77 // Build up the cell pair list. This is quite involved because we aim for a general solution that
78 // also works in narrow systems, where periodic boundary conditions make our lives hard...
79 15103 std::set<std::pair<int, int>> pairsSet; // Use set for construction to ensure uniqueness
80 // Now go through all *potential* cell pairs
81 15103 std::pair<int, int> pair;
82
2/2
✓ Branch 1 taken 407781 times.
✓ Branch 2 taken 15103 times.
 422884 for (auto &[shift1, shift2] : shiftPairs) {
83 407781 ArrayI multi1 = cellMulti + shift1;
84 407781 ArrayI multi2 = cellMulti + shift2;
85
4/4
✓ Branch 3 taken 406515 times.
✓ Branch 4 taken 1266 times.
✓ Branch 5 taken 401268 times.
✓ Branch 6 taken 5247 times.
 814296 if ((multiindex.shape == 1 && (multi1 == 1 || multi2 == 1)).any() ||
86
2/2
✓ Branch 3 taken 401268 times.
✓ Branch 4 taken 5247 times.
 406515 (multiindex.shape == 2 && (multi1 == 2 || multi2 == 2)).any()) {
87 // If there are less than 3 cells in a direction, ignore the shifts that invoke the periodic
88 // boundary conditions. This way, double counting is prevented.
89 6513 continue;
90 }
91 // Transform to linear index
92
2/2
✓ Branch 2 taken 201389 times.
✓ Branch 3 taken 199879 times.
 401268 pair = {multiindex.multi2lin(multi1, true), multiindex.multi2lin(multi2, true)};
93
3/4
✓ Branch 0 taken 201389 times.
✓ Branch 1 taken 199879 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 201389 times.
 401268 if (mode == "c08" || mode == "c04") {
94 // For newton3 modes: sort to avoid pair duplicates that only differ in order
95
2/2
✓ Branch 0 taken 96163 times.
✓ Branch 1 taken 103716 times.
 199879 auto [c1, c2] = pair;
96
2/2
✓ Branch 0 taken 96163 times.
✓ Branch 1 taken 103716 times.
 199879 if (c1 > c2) {
97 96163 pair = {c2, c1};
98 }
99 }
100
1/2
✓ Branch 1 taken 401268 times.
✗ Branch 2 not taken.
 401268 pairsSet.insert(pair);
101 }
102 // Go back to vector for random access performance
103 15103 std::vector<std::pair<int, int>> pairs;
104
2/2
✓ Branch 0 taken 303556 times.
✓ Branch 1 taken 15103 times.
 318659 for (auto &cp : pairsSet) {
105
1/2
✓ Branch 1 taken 303556 times.
✗ Branch 2 not taken.
 303556 pairs.push_back(cp);
106 }
107 30206 return pairs;
108 15103 }
109
110 108 std::vector<std::tuple<int, int, int>> Cells::getTriples(int cellIndex, const std::string &mode) {
111 108 ArrayI cellMulti = multiindex.lin2multi(cellIndex);
112
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 54 times.
 108 std::vector<std::tuple<ArrayI, ArrayI, ArrayI>> shiftTriples;
113
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 54 times.
 108 if (mode == "all") {
114 54 ArrayI shift1 = ArrayI::Zero();
115 54 ArrayI start(-1);
116 54 Multiindex<DIM> range(ArrayI(3));
117
2/2
✓ Branch 0 taken 1458 times.
✓ Branch 1 taken 54 times.
 1512 for (ArrayI shift2 : range) {
118 1458 shift2 += start;
119
2/2
✓ Branch 0 taken 39366 times.
✓ Branch 1 taken 1458 times.
 40824 for (ArrayI shift3 : range) {
120 39366 shift3 += start;
121
2/2
✓ Branch 1 taken 20844 times.
✓ Branch 2 taken 18522 times.
 39366 if (((shift2 - shift3).abs() > 1).any()) {
122 // cell 2 and 3 are not neighbors, ignore those shifts
123 20844 continue;
124 }
125
1/2
✓ Branch 1 taken 18522 times.
✗ Branch 2 not taken.
 18522 shiftTriples.push_back({shift1, shift2, shift3});
126 }
127 }
128
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
✗ Branch 3 not taken.
 54 } else if (mode == "c08" || mode == "c04") {
129 54 Multiindex<DIM> range(ArrayI(2));
130
2/2
✓ Branch 0 taken 432 times.
✓ Branch 1 taken 54 times.
 486 for (ArrayI shift1 : range) {
131
2/2
✓ Branch 0 taken 3456 times.
✓ Branch 1 taken 432 times.
 3888 for (ArrayI shift2 : range) {
132
2/2
✓ Branch 0 taken 27648 times.
✓ Branch 1 taken 3456 times.
 31104 for (ArrayI shift3 : range) {
133
2/2
✓ Branch 1 taken 9126 times.
✓ Branch 2 taken 18522 times.
 27648 if ((shift1 * shift2 * shift3 == 1).any()) {
134 9126 continue;
135 }
136
1/2
✓ Branch 1 taken 18522 times.
✗ Branch 2 not taken.
 18522 shiftTriples.push_back({shift1, shift2, shift3});
137 }
138 }
139 }
140 } else {
141 ✗ throw std::invalid_argument("Cell mode " + mode + " not implemented");
142 }
143 108 std::set<std::tuple<int, int, int>> triplesSet; // Use set for construction to ensure uniqueness
144 108 std::tuple<int, int, int> triple;
145
2/2
✓ Branch 1 taken 37044 times.
✓ Branch 2 taken 108 times.
 37152 for (auto &[shift1, shift2, shift3] : shiftTriples) {
146 37044 ArrayI multi1 = cellMulti + shift1;
147 37044 ArrayI multi2 = cellMulti + shift2;
148 37044 ArrayI multi3 = cellMulti + shift3;
149
4/4
✓ Branch 4 taken 28350 times.
✓ Branch 5 taken 8694 times.
✓ Branch 6 taken 10125 times.
✓ Branch 7 taken 18225 times.
 65394 if ((multiindex.shape == 1 && (multi1 == 1 || multi2 == 1 || multi3 == 1)).any() ||
150
2/2
✓ Branch 4 taken 10125 times.
✓ Branch 5 taken 18225 times.
 28350 (multiindex.shape == 2 && (multi1 == 2 || multi2 == 2 || multi3 == 2)).any()) {
151 // If there are less than 3 cells in a direction, ignore the shifts that invoke the periodic
152 // boundary conditions. This way, double counting is prevented.
153 26919 continue;
154 }
155 // Transform to linear index
156 10125 triple = {multiindex.multi2lin(multi1, true), multiindex.multi2lin(multi2, true),
157
2/2
✓ Branch 1 taken 4135 times.
✓ Branch 2 taken 5990 times.
 10125 multiindex.multi2lin(multi3, true)};
158 // sort to avoid duplicates that only differ in order of second and third cell
159
2/2
✓ Branch 0 taken 4135 times.
✓ Branch 1 taken 5990 times.
 10125 auto [c1, c2, c3] = triple;
160
2/2
✓ Branch 0 taken 4135 times.
✓ Branch 1 taken 5990 times.
 10125 if (c2 > c3) {
161 4135 triple = {c1, c3, c2};
162 }
163
1/2
✓ Branch 1 taken 10125 times.
✗ Branch 2 not taken.
 10125 triplesSet.insert(triple);
164 }
165 // Go back to vector for random access performance
166 108 std::vector<std::tuple<int, int, int>> triples;
167
2/2
✓ Branch 0 taken 2908 times.
✓ Branch 1 taken 108 times.
 3016 for (auto &ct : triplesSet) {
168
1/2
✓ Branch 1 taken 2908 times.
✗ Branch 2 not taken.
 2908 triples.push_back(ct);
169 }
170 216 return triples;
171 108 }
172
173 704 std::vector<int> Cells::getSameColorIndices(int color, const std::string &mode) {
174
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 704 times.
 704 std::vector<int> sameColorCells;
175
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 704 times.
✗ Branch 3 not taken.
 704 if (mode == "c08" || mode == "c04") {
176 704 std::vector<ArrayI> startingMulti;
177 704 Multiindex<DIM> rangeStartingMulti(ArrayI(2));
178
2/2
✓ Branch 0 taken 5632 times.
✓ Branch 1 taken 704 times.
 6336 for (ArrayI multi : rangeStartingMulti) {
179
1/2
✓ Branch 1 taken 5632 times.
✗ Branch 2 not taken.
 5632 startingMulti.push_back(multi);
180 }
181 704 ArrayI start = startingMulti[color];
182
2/2
✓ Branch 2 taken 130 times.
✓ Branch 3 taken 574 times.
 704 if ((multiindex.shape == 1 && start == 1).any()) {
183 // We would go around the PBC, so skip this color (i.e. leave the corresponding index list
184 // empty)
185 130 return {};
186 }
187 574 Multiindex<DIM> rangeSameColorCells(multiindex.shape - start);
188
2/2
✓ Branch 0 taken 44900 times.
✓ Branch 1 taken 574 times.
 45474 for (ArrayI multi : rangeSameColorCells) {
189
2/2
✓ Branch 0 taken 37350 times.
✓ Branch 1 taken 7550 times.
 44900 if (multi.unaryExpr([](int x) -> bool { return x % 2; }).any()) {
190 // Ignore odd differences
191 37350 continue;
192 }
193 7550 int linind = multiindex.multi2lin(start + multi, true);
194
1/2
✓ Branch 1 taken 7550 times.
✗ Branch 2 not taken.
 7550 sameColorCells.push_back(linind);
195 }
196 704 } else {
197 ✗ throw std::invalid_argument("Cell mode " + mode + " not implemented");
198 }
199 574 return sameColorCells;
200 704 }
201