• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

openmc-dev / openmc / 22914090227

10 Mar 2026 04:53PM UTC coverage: 81.569% (+0.003%) from 81.566%
22914090227

Pull #3853

github

web-flow
Merge bc1513221 into 1dc4aa988
Pull Request #3853: Fix numerical cancellation in RectLattice::distance for large pitch values

17553 of 25260 branches covered (69.49%)

Branch coverage included in aggregate %.

14 of 14 new or added lines in 2 files covered. (100.0%)

159 existing lines in 2 files now uncovered.

57964 of 67321 relevant lines covered (86.1%)

44900321.08 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

89.74
/src/lattice.cpp
1
#include "openmc/lattice.h"
2

3
#include <cmath>
4
#include <string>
5

6
#include <fmt/core.h>
7

8
#include "openmc/cell.h"
9
#include "openmc/error.h"
10
#include "openmc/geometry.h"
11
#include "openmc/geometry_aux.h"
12
#include "openmc/hdf5_interface.h"
13
#include "openmc/math_functions.h"
14
#include "openmc/string_utils.h"
15
#include "openmc/vector.h"
16
#include "openmc/xml_interface.h"
17

18
namespace openmc {
19

20
//==============================================================================
21
// Global variables
22
//==============================================================================
23

24
namespace model {
25
std::unordered_map<int32_t, int32_t> lattice_map;
26
vector<unique_ptr<Lattice>> lattices;
27
} // namespace model
28

29
//==============================================================================
30
// Lattice implementation
31
//==============================================================================
32

33
Lattice::Lattice(pugi::xml_node lat_node)
1,874✔
34
{
35
  if (check_for_node(lat_node, "id")) {
1,874!
36
    id_ = std::stoi(get_node_value(lat_node, "id"));
3,748✔
37
  } else {
38
    fatal_error("Must specify id of lattice in geometry XML file.");
×
39
  }
40

41
  if (check_for_node(lat_node, "name")) {
1,874✔
42
    name_ = get_node_value(lat_node, "name");
337✔
43
  }
44

45
  if (check_for_node(lat_node, "outer")) {
1,874✔
46
    outer_ = std::stoi(get_node_value(lat_node, "outer"));
830✔
47
  }
48
}
1,874✔
49

50
//==============================================================================
51

52
LatticeIter Lattice::begin()
289,798,651✔
53
{
54
  return LatticeIter(*this, 0);
289,798,651✔
55
}
56

57
LatticeIter Lattice::end()
48,790,595✔
58
{
59
  return LatticeIter(*this, universes_.size());
48,790,595✔
60
}
61

62
int32_t& Lattice::back()
22✔
63
{
64
  return universes_.back();
22✔
65
}
66

67
ReverseLatticeIter Lattice::rbegin()
1,638,296✔
68
{
69
  return ReverseLatticeIter(*this, universes_.size() - 1);
1,638,296✔
70
}
71

72
ReverseLatticeIter Lattice::rend()
289,047,836✔
73
{
74
  return ReverseLatticeIter(*this, -1);
289,047,836✔
75
}
76

77
//==============================================================================
78

79
void Lattice::adjust_indices()
1,874✔
80
{
81
  // Adjust the indices for the universes array.
82
  for (LatticeIter it = begin(); it != end(); ++it) {
910,194✔
83
    int uid = *it;
908,320!
84
    auto search = model::universe_map.find(uid);
908,320!
85
    if (search != model::universe_map.end()) {
908,320!
86
      *it = search->second;
908,320✔
87
    } else {
88
      fatal_error(fmt::format(
×
89
        "Invalid universe number {} specified on lattice {}", uid, id_));
×
90
    }
91
  }
92

93
  // Adjust the index for the outer universe.
94
  if (outer_ != NO_OUTER_UNIVERSE) {
1,874✔
95
    auto search = model::universe_map.find(outer_);
415!
96
    if (search != model::universe_map.end()) {
415!
97
      outer_ = search->second;
415✔
98
    } else {
99
      fatal_error(fmt::format(
×
100
        "Invalid universe number {} specified on lattice {}", outer_, id_));
×
101
    }
102
  }
103
}
1,874✔
104

105
//==============================================================================
106

107
int32_t Lattice::fill_offset_table(int32_t target_univ_id, int map,
12,609✔
108
  std::unordered_map<int32_t, int32_t>& univ_count_memo)
109
{
110
  // If the offsets have already been determined for this "map", don't bother
111
  // recalculating all of them and just return the total offset. Note that the
112
  // offsets_ array doesn't actually include the offset accounting for the last
113
  // universe, so we get the before-last offset for the given map and then
114
  // explicitly add the count for the last universe.
115
  if (offsets_[map * universes_.size() + this->begin().indx_] != C_NONE) {
12,609✔
116
    int last_offset =
33✔
117
      offsets_[(map + 1) * universes_.size() - this->begin().indx_ - 1];
33✔
118
    int last_univ = this->back();
33✔
119
    return last_offset +
33✔
120
           count_universe_instances(last_univ, target_univ_id, univ_count_memo);
33✔
121
  }
122

123
  int32_t offset = 0;
12,576✔
124
  for (LatticeIter it = begin(); it != end(); ++it) {
9,205,896✔
125
    offsets_[map * universes_.size() + it.indx_] = offset;
9,193,320✔
126
    offset += count_universe_instances(*it, target_univ_id, univ_count_memo);
9,193,320✔
127
  }
128

129
  return offset;
12,576✔
130
}
131

132
//==============================================================================
133

134
void Lattice::to_hdf5(hid_t lattices_group) const
1,432✔
135
{
136
  // Make a group for the lattice.
137
  std::string group_name {"lattice "};
1,432✔
138
  group_name += std::to_string(id_);
2,864✔
139
  hid_t lat_group = create_group(lattices_group, group_name);
1,432✔
140

141
  // Write the name and outer universe.
142
  if (!name_.empty()) {
1,432✔
143
    write_string(lat_group, "name", name_, false);
242✔
144
  }
145

146
  if (outer_ != NO_OUTER_UNIVERSE) {
1,432✔
147
    int32_t outer_id = model::universes[outer_]->id_;
316✔
148
    write_dataset(lat_group, "outer", outer_id);
316✔
149
  } else {
150
    write_dataset(lat_group, "outer", outer_);
1,116✔
151
  }
152

153
  // Call subclass-overriden function to fill in other details.
154
  to_hdf5_inner(lat_group);
1,432✔
155

156
  close_group(lat_group);
1,432✔
157
}
1,432✔
158

159
//==============================================================================
160
// RectLattice implementation
161
//==============================================================================
162

163
RectLattice::RectLattice(pugi::xml_node lat_node) : Lattice {lat_node}
1,788✔
164
{
165
  type_ = LatticeType::rect;
1,788✔
166

167
  // Read the number of lattice cells in each dimension.
168
  std::string dimension_str {get_node_value(lat_node, "dimension")};
1,788✔
169
  vector<std::string> dimension_words {split(dimension_str)};
1,788✔
170
  if (dimension_words.size() == 2) {
1,788✔
171
    n_cells_[0] = std::stoi(dimension_words[0]);
1,303✔
172
    n_cells_[1] = std::stoi(dimension_words[1]);
1,303✔
173
    n_cells_[2] = 1;
1,303✔
174
    is_3d_ = false;
1,303✔
175
  } else if (dimension_words.size() == 3) {
485!
176
    n_cells_[0] = std::stoi(dimension_words[0]);
485✔
177
    n_cells_[1] = std::stoi(dimension_words[1]);
485✔
178
    n_cells_[2] = std::stoi(dimension_words[2]);
485✔
179
    is_3d_ = true;
485✔
180
  } else {
181
    fatal_error("Rectangular lattice must be two or three dimensions.");
×
182
  }
183

184
  // Read the lattice lower-left location.
185
  std::string ll_str {get_node_value(lat_node, "lower_left")};
1,788✔
186
  vector<std::string> ll_words {split(ll_str)};
1,788✔
187
  if (ll_words.size() != dimension_words.size()) {
1,788!
188
    fatal_error("Number of entries on <lower_left> must be the same as the "
×
189
                "number of entries on <dimension>.");
190
  }
191
  lower_left_[0] = stod(ll_words[0]);
1,788✔
192
  lower_left_[1] = stod(ll_words[1]);
1,788✔
193
  if (is_3d_) {
1,788✔
194
    lower_left_[2] = stod(ll_words[2]);
485✔
195
  }
196

197
  // Read the lattice pitches.
198
  std::string pitch_str {get_node_value(lat_node, "pitch")};
1,788✔
199
  vector<std::string> pitch_words {split(pitch_str)};
1,788✔
200
  if (pitch_words.size() != dimension_words.size()) {
1,788!
201
    fatal_error("Number of entries on <pitch> must be the same as the "
×
202
                "number of entries on <dimension>.");
203
  }
204
  pitch_[0] = stod(pitch_words[0]);
1,788✔
205
  pitch_[1] = stod(pitch_words[1]);
1,788✔
206
  if (is_3d_) {
1,788✔
207
    pitch_[2] = stod(pitch_words[2]);
485✔
208
  }
209

210
  // Read the universes and make sure the correct number was specified.
211
  std::string univ_str {get_node_value(lat_node, "universes")};
1,788✔
212
  vector<std::string> univ_words {split(univ_str)};
1,788✔
213
  if (univ_words.size() != n_cells_[0] * n_cells_[1] * n_cells_[2]) {
1,788!
214
    fatal_error(fmt::format(
×
215
      "Expected {} universes for a rectangular lattice of size {}x{}x{} but {} "
216
      "were specified.",
217
      n_cells_[0] * n_cells_[1] * n_cells_[2], n_cells_[0], n_cells_[1],
×
218
      n_cells_[2], univ_words.size()));
×
219
  }
220

221
  // Parse the universes.
222
  universes_.resize(n_cells_[0] * n_cells_[1] * n_cells_[2], C_NONE);
1,788✔
223
  for (int iz = 0; iz < n_cells_[2]; iz++) {
7,967✔
224
    for (int iy = n_cells_[1] - 1; iy > -1; iy--) {
75,083✔
225
      for (int ix = 0; ix < n_cells_[0]; ix++) {
951,827✔
226
        int indx1 = n_cells_[0] * n_cells_[1] * iz +
882,923✔
227
                    n_cells_[0] * (n_cells_[1] - iy - 1) + ix;
882,923✔
228
        int indx2 = n_cells_[0] * n_cells_[1] * iz + n_cells_[0] * iy + ix;
882,923✔
229
        universes_[indx1] = std::stoi(univ_words[indx2]);
882,923✔
230
      }
231
    }
232
  }
233
}
1,788✔
234

235
//==============================================================================
236

237
const int32_t& RectLattice::operator[](const array<int, 3>& i_xyz)
859,597,662✔
238
{
239
  return universes_[get_flat_index(i_xyz)];
859,597,662✔
240
}
241

242
//==============================================================================
243

244
bool RectLattice::are_valid_indices(const array<int, 3>& i_xyz) const
2,100,265,997✔
245
{
246
  return ((i_xyz[0] >= 0) && (i_xyz[0] < n_cells_[0]) && (i_xyz[1] >= 0) &&
2,100,265,997✔
247
          (i_xyz[1] < n_cells_[1]) && (i_xyz[2] >= 0) &&
2,147,483,647!
248
          (i_xyz[2] < n_cells_[2]));
2,093,438,869!
249
}
250

251
//==============================================================================
252

253
std::pair<double, array<int, 3>> RectLattice::distance(
1,240,229,641✔
254
  Position r, Direction u, const array<int, 3>& i_xyz) const
255
{
256
  // Get short aliases to the coordinates.
257
  double x = r.x;
1,240,229,641✔
258
  double y = r.y;
1,240,229,641✔
259
  double z = r.z;
1,240,229,641✔
260

261
  // Determine the oncoming edge.
262
  double x0 {copysign(0.5 * pitch_[0], u.x)};
1,240,229,641✔
263
  double y0 {copysign(0.5 * pitch_[1], u.y)};
1,240,229,641✔
264
  double z0;
1,240,229,641✔
265

266
  // Evaluate axial distances.
267
  double dx = u.x != 0.0 ? (x0 - x) / u.x : INFTY;
1,240,229,641!
268
  double dy = u.y != 0.0 ? (y0 - y) / u.y : INFTY;
1,240,229,641!
269
  double dz = is_3d_ ? (u.z != 0.0 ? (z0 - z) / u.z : INFTY) : INFTY;
1,240,229,641!
270

271
  if (is_3d_) {
1,240,229,641✔
272
    z0 = std::copysign(0.5 * pitch_[2], u.z);
732,699,514✔
273
    dz = (u.z != 0.0) ? (z0 - z) / u.z : INFTY;
732,699,514!
274
  }
275

276
  // Minimum distance
277
  double d = std::min({dx, dy, dz});
1,240,229,641✔
278

279
  array<int, 3> lattice_trans = {0, 0, 0};
1,240,229,641✔
280

281
  // Determine which directions are crossed
282
  if (isclose(d, dx, 1e-12, FP_PRECISION))
1,240,229,641✔
283
    lattice_trans[0] = std::copysign(1, u.x);
563,894,427✔
284

285
  if (isclose(d, dy, 1e-12, FP_PRECISION))
1,240,229,641✔
286
    lattice_trans[1] = std::copysign(1, u.y);
465,946,163✔
287

288
  if (is_3d_ && isclose(d, dz, 1e-12, FP_PRECISION))
1,240,229,641✔
289
    lattice_trans[2] = std::copysign(1, u.z);
210,609,029✔
290

291
  return {d, lattice_trans};
1,240,229,641✔
292
}
293
//==============================================================================
294

295
void RectLattice::get_indices(
135,330,745✔
296
  Position r, Direction u, array<int, 3>& result) const
297
{
298
  // Determine x index, accounting for coincidence
299
  double ix_ {(r.x - lower_left_.x) / pitch_.x};
135,330,745✔
300
  long ix_close {std::lround(ix_)};
135,330,745✔
301
  if (coincident(ix_, ix_close)) {
135,330,745✔
302
    result[0] = (u.x > 0) ? ix_close : ix_close - 1;
43,893,420✔
303
  } else {
304
    result[0] = std::floor(ix_);
91,437,325✔
305
  }
306

307
  // Determine y index, accounting for coincidence
308
  double iy_ {(r.y - lower_left_.y) / pitch_.y};
135,330,745✔
309
  long iy_close {std::lround(iy_)};
135,330,745✔
310
  if (coincident(iy_, iy_close)) {
135,330,745✔
311
    result[1] = (u.y > 0) ? iy_close : iy_close - 1;
44,370,322✔
312
  } else {
313
    result[1] = std::floor(iy_);
90,960,423✔
314
  }
315

316
  // Determine z index, accounting for coincidence
317
  result[2] = 0;
135,330,745✔
318
  if (is_3d_) {
135,330,745✔
319
    double iz_ {(r.z - lower_left_.z) / pitch_.z};
69,232,847✔
320
    long iz_close {std::lround(iz_)};
69,232,847✔
321
    if (coincident(iz_, iz_close)) {
69,232,847✔
322
      result[2] = (u.z > 0) ? iz_close : iz_close - 1;
18,009,476✔
323
    } else {
324
      result[2] = std::floor(iz_);
51,223,371✔
325
    }
326
  }
327
}
135,330,745✔
328

329
int RectLattice::get_flat_index(const array<int, 3>& i_xyz) const
859,597,662✔
330
{
331
  return n_cells_[0] * n_cells_[1] * i_xyz[2] + n_cells_[0] * i_xyz[1] +
859,597,662✔
332
         i_xyz[0];
859,597,662✔
333
}
334

335
//==============================================================================
336

337
Position RectLattice::get_local_position(
863,835,071✔
338
  Position r, const array<int, 3>& i_xyz) const
339
{
340
  r.x -= (lower_left_.x + (i_xyz[0] + 0.5) * pitch_.x);
863,835,071✔
341
  r.y -= (lower_left_.y + (i_xyz[1] + 0.5) * pitch_.y);
863,835,071✔
342
  if (is_3d_) {
863,835,071✔
343
    r.z -= (lower_left_.z + (i_xyz[2] + 0.5) * pitch_.z);
683,065,862✔
344
  }
345
  return r;
863,835,071✔
346
}
347

348
//==============================================================================
349

350
int32_t& RectLattice::offset(int map, const array<int, 3>& i_xyz)
1,233,841,207✔
351
{
352
  return offsets_[n_cells_[0] * n_cells_[1] * n_cells_[2] * map +
1,233,841,207✔
353
                  n_cells_[0] * n_cells_[1] * i_xyz[2] +
1,233,841,207✔
354
                  n_cells_[0] * i_xyz[1] + i_xyz[0]];
1,233,841,207✔
355
}
356

357
//==============================================================================
358

359
int32_t RectLattice::offset(int map, int indx) const
83,295✔
360
{
361
  return offsets_[n_cells_[0] * n_cells_[1] * n_cells_[2] * map + indx];
83,295✔
362
}
363

364
//==============================================================================
365

366
std::string RectLattice::index_to_string(int indx) const
1,638,296✔
367
{
368
  int iz {indx / (n_cells_[0] * n_cells_[1])};
1,638,296✔
369
  int iy {(indx - n_cells_[0] * n_cells_[1] * iz) / n_cells_[0]};
1,638,296✔
370
  int ix {indx - n_cells_[0] * n_cells_[1] * iz - n_cells_[0] * iy};
1,638,296✔
371
  std::string out {std::to_string(ix)};
1,638,296✔
372
  out += ',';
1,638,296✔
373
  out += std::to_string(iy);
3,276,592✔
374
  if (is_3d_) {
1,638,296!
UNCOV
375
    out += ',';
×
UNCOV
376
    out += std::to_string(iz);
×
377
  }
378
  return out;
1,638,296✔
UNCOV
379
}
×
380

381
//==============================================================================
382

383
void RectLattice::to_hdf5_inner(hid_t lat_group) const
1,377✔
384
{
385
  // Write basic lattice information.
386
  write_string(lat_group, "type", "rectangular", false);
1,377✔
387
  if (is_3d_) {
1,377✔
388
    write_dataset(lat_group, "pitch", pitch_);
381✔
389
    write_dataset(lat_group, "lower_left", lower_left_);
381✔
390
    write_dataset(lat_group, "dimension", n_cells_);
381✔
391
  } else {
392
    array<double, 2> pitch_short {{pitch_[0], pitch_[1]}};
996✔
393
    write_dataset(lat_group, "pitch", pitch_short);
996✔
394
    array<double, 2> ll_short {{lower_left_[0], lower_left_[1]}};
996✔
395
    write_dataset(lat_group, "lower_left", ll_short);
996✔
396
    array<int, 2> nc_short {{n_cells_[0], n_cells_[1]}};
996✔
397
    write_dataset(lat_group, "dimension", nc_short);
996✔
398
  }
399

400
  // Write the universe ids.  The convention here is to switch the ordering on
401
  // the y-axis to match the way universes are input in a text file.
402
  if (is_3d_) {
1,377✔
403
    hsize_t nx {static_cast<hsize_t>(n_cells_[0])};
381✔
404
    hsize_t ny {static_cast<hsize_t>(n_cells_[1])};
381✔
405
    hsize_t nz {static_cast<hsize_t>(n_cells_[2])};
381✔
406
    vector<int> out(nx * ny * nz);
381✔
407

408
    for (int m = 0; m < nz; m++) {
4,061✔
409
      for (int k = 0; k < ny; k++) {
46,334✔
410
        for (int j = 0; j < nx; j++) {
568,691✔
411
          int indx1 = nx * ny * m + nx * k + j;
526,037✔
412
          int indx2 = nx * ny * m + nx * (ny - k - 1) + j;
526,037✔
413
          out[indx2] = model::universes[universes_[indx1]]->id_;
526,037✔
414
        }
415
      }
416
    }
417

418
    hsize_t dims[3] {nz, ny, nx};
381✔
419
    write_int(lat_group, 3, dims, "universes", out.data(), false);
381✔
420

421
  } else {
381✔
422
    hsize_t nx {static_cast<hsize_t>(n_cells_[0])};
996✔
423
    hsize_t ny {static_cast<hsize_t>(n_cells_[1])};
996✔
424
    vector<int> out(nx * ny);
996✔
425

426
    for (int k = 0; k < ny; k++) {
10,068✔
427
      for (int j = 0; j < nx; j++) {
149,414✔
428
        int indx1 = nx * k + j;
140,342✔
429
        int indx2 = nx * (ny - k - 1) + j;
140,342✔
430
        out[indx2] = model::universes[universes_[indx1]]->id_;
140,342✔
431
      }
432
    }
433

434
    hsize_t dims[3] {1, ny, nx};
996✔
435
    write_int(lat_group, 3, dims, "universes", out.data(), false);
996✔
436
  }
996✔
437
}
1,377✔
438

439
//==============================================================================
440
// HexLattice implementation
441
//==============================================================================
442

443
HexLattice::HexLattice(pugi::xml_node lat_node) : Lattice {lat_node}
86✔
444
{
445
  type_ = LatticeType::hex;
86✔
446

447
  // Read the number of lattice cells in each dimension.
448
  n_rings_ = std::stoi(get_node_value(lat_node, "n_rings"));
172✔
449
  if (check_for_node(lat_node, "n_axial")) {
86✔
450
    n_axial_ = std::stoi(get_node_value(lat_node, "n_axial"));
60✔
451
    is_3d_ = true;
30✔
452
  } else {
453
    n_axial_ = 1;
56✔
454
    is_3d_ = false;
56✔
455
  }
456

457
  // Read the lattice orientation.  Default to 'y'.
458
  if (check_for_node(lat_node, "orientation")) {
86✔
459
    std::string orientation = get_node_value(lat_node, "orientation");
15✔
460
    if (orientation == "y") {
15!
UNCOV
461
      orientation_ = Orientation::y;
×
462
    } else if (orientation == "x") {
15!
463
      orientation_ = Orientation::x;
15✔
464
    } else {
UNCOV
465
      fatal_error("Unrecognized orientation '" + orientation +
×
UNCOV
466
                  "' for lattice " + std::to_string(id_));
×
467
    }
468
  } else {
15✔
469
    orientation_ = Orientation::y;
71✔
470
  }
471

472
  // Read the lattice center.
473
  std::string center_str {get_node_value(lat_node, "center")};
86✔
474
  vector<std::string> center_words {split(center_str)};
86✔
475
  if (is_3d_ && (center_words.size() != 3)) {
86!
UNCOV
476
    fatal_error("A hexagonal lattice with <n_axial> must have <center> "
×
477
                "specified by 3 numbers.");
478
  } else if (!is_3d_ && center_words.size() != 2) {
86!
UNCOV
479
    fatal_error("A hexagonal lattice without <n_axial> must have <center> "
×
480
                "specified by 2 numbers.");
481
  }
482
  center_[0] = stod(center_words[0]);
86✔
483
  center_[1] = stod(center_words[1]);
86✔
484
  if (is_3d_) {
86✔
485
    center_[2] = stod(center_words[2]);
30✔
486
  }
487

488
  // Read the lattice pitches.
489
  std::string pitch_str {get_node_value(lat_node, "pitch")};
86✔
490
  vector<std::string> pitch_words {split(pitch_str)};
86✔
491
  if (is_3d_ && (pitch_words.size() != 2)) {
86!
UNCOV
492
    fatal_error("A hexagonal lattice with <n_axial> must have <pitch> "
×
493
                "specified by 2 numbers.");
494
  } else if (!is_3d_ && (pitch_words.size() != 1)) {
86!
UNCOV
495
    fatal_error("A hexagonal lattice without <n_axial> must have <pitch> "
×
496
                "specified by 1 number.");
497
  }
498
  pitch_[0] = stod(pitch_words[0]);
86✔
499
  if (is_3d_) {
86✔
500
    pitch_[1] = stod(pitch_words[1]);
30✔
501
  }
502

503
  // Read the universes and make sure the correct number was specified.
504
  int n_univ = (3 * n_rings_ * n_rings_ - 3 * n_rings_ + 1) * n_axial_;
86✔
505
  std::string univ_str {get_node_value(lat_node, "universes")};
86✔
506
  vector<std::string> univ_words {split(univ_str)};
86✔
507
  if (univ_words.size() != n_univ) {
86!
UNCOV
508
    fatal_error(fmt::format(
×
509
      "Expected {} universes for a hexagonal lattice with {} rings and {} "
510
      "axial levels but {} were specified.",
UNCOV
511
      n_univ, n_rings_, n_axial_, univ_words.size()));
×
512
  }
513

514
  // Parse the universes.
515
  // Universes in hexagonal lattices are stored in a manner that represents
516
  // a skewed coordinate system: (x, alpha) in case of 'y' orientation
517
  // and (alpha,y) in 'x' one rather than (x, y).  There is
518
  // no obvious, direct relationship between the order of universes in the
519
  // input and the order that they will be stored in the skewed array so
520
  // the following code walks a set of index values across the skewed array
521
  // in a manner that matches the input order.  Note that i_x = 0, i_a = 0
522
  // or i_a = 0, i_y = 0 corresponds to the center of the hexagonal lattice.
523
  universes_.resize((2 * n_rings_ - 1) * (2 * n_rings_ - 1) * n_axial_, C_NONE);
86✔
524
  if (orientation_ == Orientation::y) {
86✔
525
    fill_lattice_y(univ_words);
71✔
526
  } else {
527
    fill_lattice_x(univ_words);
15✔
528
  }
529
}
86✔
530

531
//==============================================================================
532

533
void HexLattice::fill_lattice_x(const vector<std::string>& univ_words)
15✔
534
{
535
  int input_index = 0;
15✔
536
  for (int m = 0; m < n_axial_; m++) {
45✔
537
    // Initialize lattice indecies.
538
    int i_a = -(n_rings_ - 1);
30✔
539
    int i_y = n_rings_ - 1;
30✔
540

541
    // Map upper region of hexagonal lattice which is found in the
542
    // first n_rings-1 rows of the input.
543
    for (int k = 0; k < n_rings_ - 1; k++) {
330✔
544

545
      // Iterate over the input columns.
546
      for (int j = 0; j < k + n_rings_; j++) {
4,950✔
547
        int indx = (2 * n_rings_ - 1) * (2 * n_rings_ - 1) * m +
4,650✔
548
                   (2 * n_rings_ - 1) * (i_y + n_rings_ - 1) +
4,650✔
549
                   (i_a + n_rings_ - 1);
4,650✔
550
        universes_[indx] = std::stoi(univ_words[input_index]);
4,650✔
551
        input_index++;
4,650✔
552
        // Move to the next right neighbour cell
553
        i_a += 1;
4,650✔
554
      }
555

556
      // Return the lattice index to the start of the current row.
557
      i_a = -(n_rings_ - 1);
300✔
558
      i_y -= 1;
300✔
559
    }
560

561
    // Map the lower region from the centerline of cart to down side
562
    for (int k = 0; k < n_rings_; k++) {
360✔
563
      // Walk the index to the lower-right neighbor of the last row start.
564
      i_a = -(n_rings_ - 1) + k;
330✔
565

566
      // Iterate over the input columns.
567
      for (int j = 0; j < 2 * n_rings_ - k - 1; j++) {
5,610✔
568
        int indx = (2 * n_rings_ - 1) * (2 * n_rings_ - 1) * m +
5,280✔
569
                   (2 * n_rings_ - 1) * (i_y + n_rings_ - 1) +
5,280✔
570
                   (i_a + n_rings_ - 1);
5,280✔
571
        universes_[indx] = std::stoi(univ_words[input_index]);
5,280✔
572
        input_index++;
5,280✔
573
        // Move to the next right neighbour cell
574
        i_a += 1;
5,280✔
575
      }
576

577
      // Return lattice index to start of current row.
578
      i_y -= 1;
330✔
579
    }
580
  }
581
}
15✔
582

583
//==============================================================================
584

585
void HexLattice::fill_lattice_y(const vector<std::string>& univ_words)
71✔
586
{
587
  int input_index = 0;
71✔
588
  for (int m = 0; m < n_axial_; m++) {
172✔
589
    // Initialize lattice indecies.
590
    int i_x = 1;
101✔
591
    int i_a = n_rings_ - 1;
101✔
592

593
    // Map upper triangular region of hexagonal lattice which is found in the
594
    // first n_rings-1 rows of the input.
595
    for (int k = 0; k < n_rings_ - 1; k++) {
622✔
596
      // Walk the index to lower-left neighbor of last row start.
597
      i_x -= 1;
521✔
598

599
      // Iterate over the input columns.
600
      for (int j = 0; j < k + 1; j++) {
3,082✔
601
        int indx = (2 * n_rings_ - 1) * (2 * n_rings_ - 1) * m +
2,561✔
602
                   (2 * n_rings_ - 1) * (i_a + n_rings_ - 1) +
2,561✔
603
                   (i_x + n_rings_ - 1);
2,561✔
604
        universes_[indx] = std::stoi(univ_words[input_index]);
2,561✔
605
        input_index++;
2,561✔
606
        // Walk the index to the right neighbor (which is not adjacent).
607
        i_x += 2;
2,561✔
608
        i_a -= 1;
2,561✔
609
      }
610

611
      // Return the lattice index to the start of the current row.
612
      i_x -= 2 * (k + 1);
521✔
613
      i_a += (k + 1);
521✔
614
    }
615

616
    // Map the middle square region of the hexagonal lattice which is found in
617
    // the next 2*n_rings-1 rows of the input.
618
    for (int k = 0; k < 2 * n_rings_ - 1; k++) {
1,244✔
619
      if ((k % 2) == 0) {
1,143✔
620
        // Walk the index to the lower-left neighbor of the last row start.
621
        i_x -= 1;
622✔
622
      } else {
623
        // Walk the index to the lower-right neighbor of the last row start.
624
        i_x += 1;
521✔
625
        i_a -= 1;
521✔
626
      }
627

628
      // Iterate over the input columns.
629
      for (int j = 0; j < n_rings_ - (k % 2); j++) {
11,488✔
630
        int indx = (2 * n_rings_ - 1) * (2 * n_rings_ - 1) * m +
10,345✔
631
                   (2 * n_rings_ - 1) * (i_a + n_rings_ - 1) +
10,345✔
632
                   (i_x + n_rings_ - 1);
10,345✔
633
        universes_[indx] = std::stoi(univ_words[input_index]);
10,345✔
634
        input_index++;
10,345✔
635
        // Walk the index to the right neighbor (which is not adjacent).
636
        i_x += 2;
10,345✔
637
        i_a -= 1;
10,345✔
638
      }
639

640
      // Return the lattice index to the start of the current row.
641
      i_x -= 2 * (n_rings_ - (k % 2));
1,143✔
642
      i_a += n_rings_ - (k % 2);
1,143✔
643
    }
644

645
    // Map the lower triangular region of the hexagonal lattice.
646
    for (int k = 0; k < n_rings_ - 1; k++) {
622✔
647
      // Walk the index to the lower-right neighbor of the last row start.
648
      i_x += 1;
521✔
649
      i_a -= 1;
521✔
650

651
      // Iterate over the input columns.
652
      for (int j = 0; j < n_rings_ - k - 1; j++) {
3,082✔
653
        int indx = (2 * n_rings_ - 1) * (2 * n_rings_ - 1) * m +
2,561✔
654
                   (2 * n_rings_ - 1) * (i_a + n_rings_ - 1) +
2,561✔
655
                   (i_x + n_rings_ - 1);
2,561✔
656
        universes_[indx] = std::stoi(univ_words[input_index]);
2,561✔
657
        input_index++;
2,561✔
658
        // Walk the index to the right neighbor (which is not adjacent).
659
        i_x += 2;
2,561✔
660
        i_a -= 1;
2,561✔
661
      }
662

663
      // Return lattice index to start of current row.
664
      i_x -= 2 * (n_rings_ - k - 1);
521✔
665
      i_a += n_rings_ - k - 1;
521✔
666
    }
667
  }
668
}
71✔
669

670
//==============================================================================
671

672
const int32_t& HexLattice::operator[](const array<int, 3>& i_xyz)
29,285,960✔
673
{
674
  return universes_[get_flat_index(i_xyz)];
29,285,960✔
675
}
676

677
//==============================================================================
678

679
// The HexLattice iterators need their own versions b/c the universes array is
680
// "square", meaning that it is allocated with entries that are intentionally
681
// left empty. As such, the iterator indices need to skip the empty entries to
682
// get cell instances and geometry paths correct. See the image in the Theory
683
// and Methodology section on "Hexagonal Lattice Indexing" for a visual of where
684
// the empty positions are.
685
LatticeIter HexLattice::begin()
2,086✔
686
{
687
  return LatticeIter(*this, n_rings_ - 1);
2,086✔
688
}
689

690
ReverseLatticeIter HexLattice::rbegin()
231✔
691
{
692
  return ReverseLatticeIter(*this, universes_.size() - n_rings_);
231✔
693
}
694

695
int32_t& HexLattice::back()
11✔
696
{
697
  return universes_[universes_.size() - n_rings_];
11✔
698
}
699

700
LatticeIter HexLattice::end()
643,129✔
701
{
702
  return LatticeIter(*this, universes_.size() - n_rings_ + 1);
643,129✔
703
}
704

705
ReverseLatticeIter HexLattice::rend()
924✔
706
{
707
  return ReverseLatticeIter(*this, n_rings_ - 2);
924✔
708
}
709

710
//==============================================================================
711

712
bool HexLattice::are_valid_indices(const array<int, 3>& i_xyz) const
113,638,691✔
713
{
714
  // Check if (x, alpha, z) indices are valid, accounting for number of rings
715
  return ((i_xyz[0] >= 0) && (i_xyz[1] >= 0) && (i_xyz[2] >= 0) &&
113,638,691✔
716
          (i_xyz[0] < 2 * n_rings_ - 1) && (i_xyz[1] < 2 * n_rings_ - 1) &&
110,318,660✔
717
          (i_xyz[0] + i_xyz[1] > n_rings_ - 2) &&
109,545,008✔
718
          (i_xyz[0] + i_xyz[1] < 3 * n_rings_ - 2) && (i_xyz[2] < n_axial_));
222,837,269✔
719
}
720

721
//==============================================================================
722

723
std::pair<double, array<int, 3>> HexLattice::distance(
85,503,297✔
724
  Position r, Direction u, const array<int, 3>& i_xyz) const
725
{
726
  // Short description of the direction vectors used here.  The beta, gamma, and
727
  // delta vectors point towards the flat sides of each hexagonal tile.
728
  // Y - orientation:
729
  //   basis0 = (1, 0)
730
  //   basis1 = (-1/sqrt(3), 1)   = +120 degrees from basis0
731
  //   beta   = (sqrt(3)/2, 1/2)  = +30 degrees from basis0
732
  //   gamma  = (sqrt(3)/2, -1/2) = -60 degrees from beta
733
  //   delta  = (0, 1)            = +60 degrees from beta
734
  // X - orientation:
735
  //   basis0 = (1/sqrt(3), -1)
736
  //   basis1 = (0, 1)            = +120 degrees from basis0
737
  //   beta   = (1, 0)            = +30 degrees from basis0
738
  //   gamma  = (1/2, -sqrt(3)/2) = -60 degrees from beta
739
  //   delta  = (1/2, sqrt(3)/2)  = +60 degrees from beta
740
  // The z-axis is considered separately.
741
  double beta_dir;
85,503,297✔
742
  double gamma_dir;
85,503,297✔
743
  double delta_dir;
85,503,297✔
744
  if (orientation_ == Orientation::y) {
85,503,297✔
745
    beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0;
70,869,348✔
746
    gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0;
70,869,348✔
747
    delta_dir = u.y;
70,869,348✔
748
  } else {
749
    beta_dir = u.x;
14,633,949✔
750
    gamma_dir = u.x / 2.0 - u.y * std::sqrt(3.0) / 2.0;
14,633,949✔
751
    delta_dir = u.x / 2.0 + u.y * std::sqrt(3.0) / 2.0;
14,633,949✔
752
  }
753

754
  // Note that hexagonal lattice distance calculations are performed
755
  // using the particle's coordinates relative to the neighbor lattice
756
  // cells, not relative to the particle's current cell.  This is done
757
  // because there is significant disagreement between neighboring cells
758
  // on where the lattice boundary is due to finite precision issues.
759

760
  // beta direction
761
  double d {INFTY};
85,503,297✔
762
  array<int, 3> lattice_trans;
85,503,297✔
763
  double edge = -copysign(0.5 * pitch_[0], beta_dir); // Oncoming edge
85,503,297✔
764
  Position r_t;
85,503,297✔
765
  if (beta_dir > 0) {
85,503,297✔
766
    const array<int, 3> i_xyz_t {i_xyz[0] + 1, i_xyz[1], i_xyz[2]};
42,766,713✔
767
    r_t = get_local_position(r, i_xyz_t);
42,766,713✔
768
  } else {
769
    const array<int, 3> i_xyz_t {i_xyz[0] - 1, i_xyz[1], i_xyz[2]};
42,736,584✔
770
    r_t = get_local_position(r, i_xyz_t);
42,736,584✔
771
  }
772
  double beta;
85,503,297✔
773
  if (orientation_ == Orientation::y) {
85,503,297✔
774
    beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0;
70,869,348✔
775
  } else {
776
    beta = r_t.x;
777
  }
778
  if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) {
85,503,297!
779
    d = (edge - beta) / beta_dir;
85,503,297✔
780
    if (beta_dir > 0) {
85,503,297✔
781
      lattice_trans = {1, 0, 0};
782
    } else {
783
      lattice_trans = {-1, 0, 0};
42,736,584✔
784
    }
785
  }
786

787
  // gamma direction
788
  edge = -copysign(0.5 * pitch_[0], gamma_dir);
85,503,297✔
789
  if (gamma_dir > 0) {
85,503,297✔
790
    const array<int, 3> i_xyz_t {i_xyz[0] + 1, i_xyz[1] - 1, i_xyz[2]};
42,833,450✔
791
    r_t = get_local_position(r, i_xyz_t);
42,833,450✔
792
  } else {
793
    const array<int, 3> i_xyz_t {i_xyz[0] - 1, i_xyz[1] + 1, i_xyz[2]};
42,669,847✔
794
    r_t = get_local_position(r, i_xyz_t);
42,669,847✔
795
  }
796
  double gamma;
85,503,297✔
797
  if (orientation_ == Orientation::y) {
85,503,297✔
798
    gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0;
70,869,348✔
799
  } else {
800
    gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0;
14,633,949✔
801
  }
802
  if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) {
85,503,297!
803
    double this_d = (edge - gamma) / gamma_dir;
85,503,297✔
804
    if (this_d < d) {
85,503,297✔
805
      if (gamma_dir > 0) {
42,791,749✔
806
        lattice_trans = {1, -1, 0};
807
      } else {
808
        lattice_trans = {-1, 1, 0};
21,307,407✔
809
      }
810
      d = this_d;
811
    }
812
  }
813

814
  // delta direction
815
  edge = -copysign(0.5 * pitch_[0], delta_dir);
85,503,297✔
816
  if (delta_dir > 0) {
85,503,297✔
817
    const array<int, 3> i_xyz_t {i_xyz[0], i_xyz[1] + 1, i_xyz[2]};
42,674,093✔
818
    r_t = get_local_position(r, i_xyz_t);
42,674,093✔
819
  } else {
820
    const array<int, 3> i_xyz_t {i_xyz[0], i_xyz[1] - 1, i_xyz[2]};
42,829,204✔
821
    r_t = get_local_position(r, i_xyz_t);
42,829,204✔
822
  }
823
  double delta;
85,503,297✔
824
  if (orientation_ == Orientation::y) {
85,503,297✔
825
    delta = r_t.y;
826
  } else {
827
    delta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0;
14,633,949✔
828
  }
829
  if ((std::abs(delta - edge) > FP_PRECISION) && delta_dir != 0) {
85,503,297!
830
    double this_d = (edge - delta) / delta_dir;
85,503,297✔
831
    if (this_d < d) {
85,503,297✔
832
      if (delta_dir > 0) {
28,519,040✔
833
        lattice_trans = {0, 1, 0};
834
      } else {
835
        lattice_trans = {0, -1, 0};
14,284,358✔
836
      }
837
      d = this_d;
838
    }
839
  }
840

841
  // Top and bottom sides
842
  if (is_3d_) {
85,503,297✔
843
    double z = r.z;
20,213,028✔
844
    double z0 {copysign(0.5 * pitch_[1], u.z)};
20,213,028!
845
    if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) {
20,213,028!
846
      double this_d = (z0 - z) / u.z;
20,213,028✔
847
      if (this_d < d) {
20,213,028✔
848
        d = this_d;
3,062,268✔
849
        if (u.z > 0) {
3,062,268✔
850
          lattice_trans = {0, 0, 1};
851
        } else {
852
          lattice_trans = {0, 0, -1};
1,535,215✔
853
        }
854
        d = this_d;
855
      }
856
    }
857
  }
858

859
  return {d, lattice_trans};
85,503,297✔
860
}
861

862
//==============================================================================
863

864
void HexLattice::get_indices(
11,900,471✔
865
  Position r, Direction u, array<int, 3>& result) const
866
{
867
  // Offset the xyz by the lattice center.
868
  Position r_o {r.x - center_.x, r.y - center_.y, r.z};
11,900,471✔
869
  if (is_3d_) {
11,900,471✔
870
    r_o.z -= center_.z;
2,484,878✔
871
  }
872

873
  // Index the z direction, accounting for coincidence
874
  result[2] = 0;
11,900,471✔
875
  if (is_3d_) {
11,900,471✔
876
    double iz_ {r_o.z / pitch_[1] + 0.5 * n_axial_};
2,484,878✔
877
    long iz_close {std::lround(iz_)};
2,484,878✔
878
    if (coincident(iz_, iz_close)) {
2,484,878✔
879
      result[2] = (u.z > 0) ? iz_close : iz_close - 1;
765,699✔
880
    } else {
881
      result[2] = std::floor(iz_);
1,719,179✔
882
    }
883
  }
884

885
  if (orientation_ == Orientation::y) {
11,900,471✔
886
    // Convert coordinates into skewed bases.  The (x, alpha) basis is used to
887
    // find the index of the global coordinates to within 4 cells.
888
    double alpha = r_o.y - r_o.x / std::sqrt(3.0);
11,078,199✔
889
    result[0] = std::floor(r_o.x / (0.5 * std::sqrt(3.0) * pitch_[0]));
11,078,199✔
890
    result[1] = std::floor(alpha / pitch_[0]);
11,078,199✔
891
  } else {
892
    // Convert coordinates into skewed bases.  The (alpha, y) basis is used to
893
    // find the index of the global coordinates to within 4 cells.
894
    double alpha = r_o.y - r_o.x * std::sqrt(3.0);
822,272✔
895
    result[0] = std::floor(-alpha / (std::sqrt(3.0) * pitch_[0]));
822,272✔
896
    result[1] = std::floor(r_o.y / (0.5 * std::sqrt(3.0) * pitch_[0]));
822,272✔
897
  }
898

899
  // Add offset to indices (the center cell is (i1, i2) = (0, 0) but
900
  // the array is offset so that the indices never go below 0).
901
  result[0] += n_rings_ - 1;
11,900,471✔
902
  result[1] += n_rings_ - 1;
11,900,471✔
903

904
  // Calculate the (squared) distance between the particle and the centers of
905
  // the four possible cells.  Regular hexagonal tiles form a Voronoi
906
  // tessellation so the xyz should be in the hexagonal cell that it is closest
907
  // to the center of.  This method is used over a method that uses the
908
  // remainders of the floor divisions above because it provides better finite
909
  // precision performance.  Squared distances are used because they are more
910
  // computationally efficient than normal distances.
911

912
  // COINCIDENCE CHECK
913
  // if a distance to center, d, is within the coincidence tolerance of the
914
  // current minimum distance, d_min, the particle is on an edge or vertex.
915
  // In this case, the dot product of the position vector and direction vector
916
  // for the current indices, dp, and the dot product for the currently selected
917
  // indices, dp_min, are compared. The cell which the particle is moving into
918
  // is kept (i.e. the cell with the lowest dot product as the vectors will be
919
  // completely opposed if the particle is moving directly toward the center of
920
  // the cell).
921
  int i1_chg;
11,900,471✔
922
  int i2_chg;
11,900,471✔
923
  double d_min {INFTY};
11,900,471✔
924
  double dp_min {INFTY};
11,900,471✔
925
  for (int i = 0; i < 2; i++) {
35,701,413✔
926
    for (int j = 0; j < 2; j++) {
71,402,826✔
927
      // get local coordinates
928
      const array<int, 3> i_xyz {result[0] + j, result[1] + i, 0};
47,601,884✔
929
      Position r_t = get_local_position(r, i_xyz);
47,601,884✔
930
      // calculate distance
931
      double d = r_t.x * r_t.x + r_t.y * r_t.y;
47,601,884✔
932
      // check for coincidence. Because the numerical error incurred
933
      // in hex geometry is higher than other geometries, the relative
934
      // coincidence is checked here so that coincidence is successfully
935
      // detected on large hex lattice with particles far from the origin
936
      // which have rounding errors larger than the FP_COINCIDENT thresdhold.
937
      bool on_boundary = coincident(1.0, d_min / d);
47,601,884✔
938
      if (d < d_min || on_boundary) {
47,601,884✔
939
        // normalize r_t and find dot product
940
        r_t /= std::sqrt(d);
29,822,133✔
941
        double dp = u.x * r_t.x + u.y * r_t.y;
29,822,133✔
942
        // do not update values if particle is on a
943
        // boundary and not moving into this cell
944
        if (on_boundary && dp > dp_min)
29,822,133✔
945
          continue;
2,050,257✔
946
        // update values
947
        d_min = d;
948
        i1_chg = j;
949
        i2_chg = i;
950
        dp_min = dp;
951
      }
952
    }
953
  }
954

955
  // update outgoing indices
956
  result[0] += i1_chg;
11,900,471✔
957
  result[1] += i2_chg;
11,900,471✔
958
}
11,900,471✔
959

960
int HexLattice::get_flat_index(const array<int, 3>& i_xyz) const
29,285,960✔
961
{
962
  return (2 * n_rings_ - 1) * (2 * n_rings_ - 1) * i_xyz[2] +
29,285,960✔
963
         (2 * n_rings_ - 1) * i_xyz[1] + i_xyz[0];
29,285,960✔
964
}
965

966
//==============================================================================
967

968
Position HexLattice::get_local_position(
337,758,718✔
969
  Position r, const array<int, 3>& i_xyz) const
970
{
971
  if (orientation_ == Orientation::y) {
337,758,718✔
972
    // x_l = x_g - (center + pitch_x*cos(30)*index_x)
973
    r.x -=
572,539,000✔
974
      center_.x + std::sqrt(3.0) / 2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0];
286,269,500✔
975
    // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y)
976
    r.y -= (center_.y + (i_xyz[1] - n_rings_ + 1) * pitch_[0] +
286,269,500✔
977
            (i_xyz[0] - n_rings_ + 1) * pitch_[0] / 2.0);
286,269,500✔
978
  } else {
979
    // x_l = x_g - (center + pitch_x*index_a + pitch_y*sin(30)*index_y)
980
    r.x -= (center_.x + (i_xyz[0] - n_rings_ + 1) * pitch_[0] +
51,489,218✔
981
            (i_xyz[1] - n_rings_ + 1) * pitch_[0] / 2.0);
51,489,218✔
982
    // y_l = y_g - (center + pitch_y*cos(30)*index_y)
983
    r.y -=
51,489,218✔
984
      center_.y + std::sqrt(3.0) / 2.0 * (i_xyz[1] - n_rings_ + 1) * pitch_[0];
51,489,218✔
985
  }
986

987
  if (is_3d_) {
337,758,718✔
988
    r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1];
78,611,665✔
989
  }
990

991
  return r;
337,758,718✔
992
}
993

994
//==============================================================================
995

996
bool HexLattice::is_valid_index(int indx) const
363,287✔
997
{
998
  int nx {2 * n_rings_ - 1};
363,287✔
999
  int ny {2 * n_rings_ - 1};
363,287✔
1000
  int iz = indx / (nx * ny);
363,287✔
1001
  int iy = (indx - nx * ny * iz) / nx;
363,287✔
1002
  int ix = indx - nx * ny * iz - nx * iy;
363,287✔
1003
  array<int, 3> i_xyz {ix, iy, iz};
363,287✔
1004
  return are_valid_indices(i_xyz);
363,287✔
1005
}
1006

1007
//==============================================================================
1008

1009
int32_t& HexLattice::offset(int map, const array<int, 3>& i_xyz)
77,204,831✔
1010
{
1011
  int nx {2 * n_rings_ - 1};
77,204,831✔
1012
  int ny {2 * n_rings_ - 1};
77,204,831✔
1013
  int nz {n_axial_};
77,204,831✔
1014
  return offsets_[nx * ny * nz * map + nx * ny * i_xyz[2] + nx * i_xyz[1] +
77,204,831✔
1015
                  i_xyz[0]];
77,204,831✔
1016
}
1017

UNCOV
1018
int32_t HexLattice::offset(int map, int indx) const
×
1019
{
1020
  int nx {2 * n_rings_ - 1};
×
UNCOV
1021
  int ny {2 * n_rings_ - 1};
×
1022
  int nz {n_axial_};
×
1023
  return offsets_[nx * ny * nz * map + indx];
×
1024
}
1025

1026
//==============================================================================
1027

1028
std::string HexLattice::index_to_string(int indx) const
231✔
1029
{
1030
  int nx {2 * n_rings_ - 1};
231✔
1031
  int ny {2 * n_rings_ - 1};
231✔
1032
  int iz {indx / (nx * ny)};
231✔
1033
  int iy {(indx - nx * ny * iz) / nx};
231✔
1034
  int ix {indx - nx * ny * iz - nx * iy};
231✔
1035
  std::string out {std::to_string(ix - n_rings_ + 1)};
231✔
1036
  out += ',';
231✔
1037
  out += std::to_string(iy - n_rings_ + 1);
462✔
1038
  if (is_3d_) {
231!
UNCOV
1039
    out += ',';
×
UNCOV
1040
    out += std::to_string(iz);
×
1041
  }
1042
  return out;
231✔
UNCOV
1043
}
×
1044

1045
//==============================================================================
1046

1047
void HexLattice::to_hdf5_inner(hid_t lat_group) const
55✔
1048
{
1049
  // Write basic lattice information.
1050
  write_string(lat_group, "type", "hexagonal", false);
55✔
1051
  write_dataset(lat_group, "n_rings", n_rings_);
55✔
1052
  write_dataset(lat_group, "n_axial", n_axial_);
55✔
1053
  if (orientation_ == Orientation::y) {
55✔
1054
    write_string(lat_group, "orientation", "y", false);
88✔
1055
  } else {
1056
    write_string(lat_group, "orientation", "x", false);
22✔
1057
  }
1058
  if (is_3d_) {
55✔
1059
    write_dataset(lat_group, "pitch", pitch_);
22✔
1060
    write_dataset(lat_group, "center", center_);
22✔
1061
  } else {
1062
    array<double, 1> pitch_short {{pitch_[0]}};
33✔
1063
    write_dataset(lat_group, "pitch", pitch_short);
33✔
1064
    array<double, 2> center_short {{center_[0], center_[1]}};
33✔
1065
    write_dataset(lat_group, "center", center_short);
33✔
1066
  }
1067

1068
  // Write the universe ids.
1069
  hsize_t nx {static_cast<hsize_t>(2 * n_rings_ - 1)};
55✔
1070
  hsize_t ny {static_cast<hsize_t>(2 * n_rings_ - 1)};
55✔
1071
  hsize_t nz {static_cast<hsize_t>(n_axial_)};
55✔
1072
  vector<int> out(nx * ny * nz);
55✔
1073

1074
  for (int m = 0; m < nz; m++) {
143✔
1075
    for (int k = 0; k < ny; k++) {
1,364✔
1076
      for (int j = 0; j < nx; j++) {
26,004✔
1077
        int indx = nx * ny * m + nx * k + j;
24,728✔
1078
        if (j + k < n_rings_ - 1) {
24,728✔
1079
          // This array position is never used; put a -1 to indicate this.
1080
          out[indx] = -1;
3,080✔
1081
        } else if (j + k > 3 * n_rings_ - 3) {
21,648✔
1082
          // This array position is never used; put a -1 to indicate this.
1083
          out[indx] = -1;
3,080✔
1084
        } else {
1085
          out[indx] = model::universes[universes_[indx]]->id_;
18,568✔
1086
        }
1087
      }
1088
    }
1089
  }
1090

1091
  hsize_t dims[3] {nz, ny, nx};
55✔
1092
  write_int(lat_group, 3, dims, "universes", out.data(), false);
55✔
1093
}
55✔
1094

1095
//==============================================================================
1096
// Non-method functions
1097
//==============================================================================
1098

1099
void read_lattices(pugi::xml_node node)
8,156✔
1100
{
1101
  for (pugi::xml_node lat_node : node.children("lattice")) {
9,944✔
1102
    model::lattices.push_back(make_unique<RectLattice>(lat_node));
1,788✔
1103
  }
1104
  for (pugi::xml_node lat_node : node.children("hex_lattice")) {
8,242✔
1105
    model::lattices.push_back(make_unique<HexLattice>(lat_node));
86✔
1106
  }
1107

1108
  // Fill the lattice map.
1109
  for (int i_lat = 0; i_lat < model::lattices.size(); i_lat++) {
10,030✔
1110
    int id = model::lattices[i_lat]->id_;
1,874!
1111
    auto in_map = model::lattice_map.find(id);
1,874!
1112
    if (in_map == model::lattice_map.end()) {
1,874!
1113
      model::lattice_map[id] = i_lat;
1,874✔
1114
    } else {
UNCOV
1115
      fatal_error(
×
UNCOV
1116
        fmt::format("Two or more lattices use the same unique ID: {}", id));
×
1117
    }
1118
  }
1119
}
8,156✔
1120

1121
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc