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 |