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

openmc-dev / openmc / 27291283994

10 Jun 2026 04:44PM UTC coverage: 81.276% (-0.005%) from 81.281%
27291283994

Pull #3934

github

web-flow
Merge 4a81ca18b into 02eb999af
Pull Request #3934: Fix virtual surface crossing

18153 of 26343 branches covered (68.91%)

Branch coverage included in aggregate %.

27 of 27 new or added lines in 1 file covered. (100.0%)

1 existing line in 1 file now uncovered.

59271 of 68918 relevant lines covered (86.0%)

48701050.65 hits per line

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

75.14
/src/cell.cpp
1

2
#include "openmc/cell.h"
3

4
#include <algorithm>
5
#include <cassert>
6
#include <cctype>
7
#include <cmath>
8
#include <iterator>
9
#include <set>
10
#include <sstream>
11
#include <string>
12

13
#include <fmt/core.h>
14

15
#include "openmc/capi.h"
16
#include "openmc/constants.h"
17
#include "openmc/dagmc.h"
18
#include "openmc/error.h"
19
#include "openmc/geometry.h"
20
#include "openmc/hdf5_interface.h"
21
#include "openmc/lattice.h"
22
#include "openmc/material.h"
23
#include "openmc/nuclide.h"
24
#include "openmc/settings.h"
25
#include "openmc/xml_interface.h"
26

27
namespace openmc {
28

29
//==============================================================================
30
// Global variables
31
//==============================================================================
32

33
namespace model {
34
std::unordered_map<int32_t, int32_t> cell_map;
35
vector<unique_ptr<Cell>> cells;
36

37
} // namespace model
38

39
//==============================================================================
40
// Cell implementation
41
//==============================================================================
42

43
int32_t Cell::n_instances() const
15,760✔
44
{
45
  return model::universes[universe_]->n_instances_;
15,760✔
46
}
47

48
void Cell::set_rotation(const vector<double>& rot)
445✔
49
{
50
  if (fill_ == C_NONE) {
445!
51
    fatal_error(fmt::format("Cannot apply a rotation to cell {}"
×
52
                            " because it is not filled with another universe",
53
      id_));
×
54
  }
55

56
  if (rot.size() != 3 && rot.size() != 9) {
445!
57
    fatal_error(fmt::format("Non-3D rotation vector applied to cell {}", id_));
×
58
  }
59

60
  // Compute and store the inverse rotation matrix for the angles given.
61
  rotation_.clear();
445✔
62
  rotation_.reserve(rot.size() == 9 ? 9 : 12);
890!
63
  if (rot.size() == 3) {
445!
64
    double phi = -rot[0] * PI / 180.0;
445✔
65
    double theta = -rot[1] * PI / 180.0;
445✔
66
    double psi = -rot[2] * PI / 180.0;
445✔
67
    rotation_.push_back(std::cos(theta) * std::cos(psi));
445✔
68
    rotation_.push_back(-std::cos(phi) * std::sin(psi) +
445✔
69
                        std::sin(phi) * std::sin(theta) * std::cos(psi));
445✔
70
    rotation_.push_back(std::sin(phi) * std::sin(psi) +
445✔
71
                        std::cos(phi) * std::sin(theta) * std::cos(psi));
445✔
72
    rotation_.push_back(std::cos(theta) * std::sin(psi));
445✔
73
    rotation_.push_back(std::cos(phi) * std::cos(psi) +
445✔
74
                        std::sin(phi) * std::sin(theta) * std::sin(psi));
445✔
75
    rotation_.push_back(-std::sin(phi) * std::cos(psi) +
445✔
76
                        std::cos(phi) * std::sin(theta) * std::sin(psi));
445✔
77
    rotation_.push_back(-std::sin(theta));
445✔
78
    rotation_.push_back(std::sin(phi) * std::cos(theta));
445✔
79
    rotation_.push_back(std::cos(phi) * std::cos(theta));
445✔
80

81
    // When user specifies angles, write them at end of vector
82
    rotation_.push_back(rot[0]);
445✔
83
    rotation_.push_back(rot[1]);
445✔
84
    rotation_.push_back(rot[2]);
445✔
85
  } else {
86
    std::copy(rot.begin(), rot.end(), std::back_inserter(rotation_));
×
87
  }
88
}
445✔
89

90
double Cell::temperature(int32_t instance) const
9,634✔
91
{
92
  if (sqrtkT_.size() < 1) {
9,634!
93
    throw std::runtime_error {"Cell temperature has not yet been set."};
×
94
  }
95

96
  if (instance >= 0) {
9,634✔
97
    double sqrtkT = sqrtkT_.size() == 1 ? sqrtkT_.at(0) : sqrtkT_.at(instance);
9,548✔
98
    return sqrtkT * sqrtkT / K_BOLTZMANN;
9,548✔
99
  } else {
100
    return sqrtkT_[0] * sqrtkT_[0] / K_BOLTZMANN;
86✔
101
  }
102
}
103

104
double Cell::density_mult(int32_t instance) const
2,147,483,647✔
105
{
106
  if (instance >= 0) {
2,147,483,647✔
107
    return density_mult_.size() == 1 ? density_mult_.at(0)
2,147,483,647✔
108
                                     : density_mult_.at(instance);
5,034,975✔
109
  } else {
110
    return density_mult_[0];
77✔
111
  }
112
}
113

114
double Cell::density(int32_t instance) const
198✔
115
{
116
  const int32_t mat_index = material(instance);
198!
117
  if (mat_index == MATERIAL_VOID)
198!
118
    return 0.0;
119

120
  return density_mult(instance) * model::materials[mat_index]->density_gpcc();
396!
121
}
122

123
void Cell::set_temperature(double T, int32_t instance, bool set_contained)
10,012✔
124
{
125
  if (settings::temperature_method == TemperatureMethod::INTERPOLATION) {
10,012!
126
    if (T < (data::temperature_min - settings::temperature_tolerance)) {
×
127
      throw std::runtime_error {
×
128
        fmt::format("Temperature of {} K is below minimum temperature at "
×
129
                    "which data is available of {} K.",
130
          T, data::temperature_min)};
×
131
    } else if (T > (data::temperature_max + settings::temperature_tolerance)) {
×
132
      throw std::runtime_error {
×
133
        fmt::format("Temperature of {} K is above maximum temperature at "
×
134
                    "which data is available of {} K.",
135
          T, data::temperature_max)};
×
136
    }
137
  }
138

139
  if (type_ == Fill::MATERIAL) {
10,012✔
140
    if (instance >= 0) {
9,982✔
141
      // If temperature vector is not big enough, resize it first
142
      if (sqrtkT_.size() != n_instances())
9,905✔
143
        sqrtkT_.resize(n_instances(), sqrtkT_[0]);
45✔
144

145
      // Set temperature for the corresponding instance
146
      sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T);
9,905✔
147
    } else {
148
      // Set temperature for all instances
149
      for (auto& T_ : sqrtkT_) {
154✔
150
        T_ = std::sqrt(K_BOLTZMANN * T);
77✔
151
      }
152
    }
153
  } else {
154
    if (!set_contained) {
30!
155
      throw std::runtime_error {
×
156
        fmt::format("Attempted to set the temperature of cell {} "
×
157
                    "which is not filled by a material.",
158
          id_)};
×
159
    }
160

161
    auto contained_cells = this->get_contained_cells(instance);
30✔
162
    for (const auto& entry : contained_cells) {
120✔
163
      auto& cell = model::cells[entry.first];
90!
164
      assert(cell->type_ == Fill::MATERIAL);
90!
165
      auto& instances = entry.second;
90✔
166
      for (auto instance : instances) {
315✔
167
        cell->set_temperature(T, instance);
225✔
168
      }
169
    }
170
  }
30✔
171
}
10,012✔
172

173
void Cell::set_density(double density, int32_t instance, bool set_contained)
346✔
174
{
175
  if (type_ != Fill::MATERIAL && !set_contained) {
346!
176
    fatal_error(
×
177
      fmt::format("Attempted to set the density multiplier of cell {} "
×
178
                  "which is not filled by a material.",
179
        id_));
×
180
  }
181

182
  if (type_ == Fill::MATERIAL) {
346✔
183
    const int32_t mat_index = material(instance);
331!
184
    if (mat_index == MATERIAL_VOID)
331!
185
      return;
186

187
    if (instance >= 0) {
331✔
188
      // If density multiplier vector is not big enough, resize it first
189
      if (density_mult_.size() != n_instances())
254✔
190
        density_mult_.resize(n_instances(), density_mult_[0]);
111✔
191

192
      // Set density multiplier for the corresponding instance
193
      density_mult_.at(instance) =
254✔
194
        density / model::materials[mat_index]->density_gpcc();
508!
195
    } else {
196
      // Set density multiplier for all instances
197
      for (auto& x : density_mult_) {
154✔
198
        x = density / model::materials[mat_index]->density_gpcc();
154!
199
      }
200
    }
201
  } else {
202
    auto contained_cells = this->get_contained_cells(instance);
15✔
203
    for (const auto& entry : contained_cells) {
60✔
204
      auto& cell = model::cells[entry.first];
45!
205
      assert(cell->type_ == Fill::MATERIAL);
45!
206
      auto& instances = entry.second;
45✔
207
      for (auto instance : instances) {
90✔
208
        cell->set_density(density, instance);
45✔
209
      }
210
    }
211
  }
15✔
212
}
213

214
void Cell::export_properties_hdf5(hid_t group) const
231✔
215
{
216
  // Create a group for this cell.
217
  auto cell_group = create_group(group, fmt::format("cell {}", id_));
231✔
218

219
  // Write temperature in [K] for one or more cell instances
220
  vector<double> temps;
231✔
221
  for (auto sqrtkT_val : sqrtkT_)
429✔
222
    temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN);
198✔
223
  write_dataset(cell_group, "temperature", temps);
231✔
224

225
  // Write density for one or more cell instances
226
  if (type_ == Fill::MATERIAL && material_.size() > 0) {
231✔
227
    vector<double> density;
198✔
228
    for (int32_t i = 0; i < density_mult_.size(); ++i)
396✔
229
      density.push_back(this->density(i));
198✔
230

231
    write_dataset(cell_group, "density", density);
198✔
232
  }
198✔
233

234
  close_group(cell_group);
231✔
235
}
231✔
236

237
void Cell::import_properties_hdf5(hid_t group)
253✔
238
{
239
  auto cell_group = open_group(group, fmt::format("cell {}", id_));
253✔
240

241
  // Read temperatures from file
242
  vector<double> temps;
253✔
243
  read_dataset(cell_group, "temperature", temps);
253✔
244

245
  // Ensure number of temperatures makes sense
246
  auto n_temps = temps.size();
253✔
247
  if (n_temps > 1 && n_temps != n_instances()) {
253!
248
    fatal_error(fmt::format(
×
249
      "Number of temperatures for cell {} doesn't match number of instances",
250
      id_));
×
251
  }
252

253
  // Modify temperatures for the cell
254
  sqrtkT_.clear();
253✔
255
  sqrtkT_.resize(temps.size());
253✔
256
  for (int64_t i = 0; i < temps.size(); ++i) {
9,922✔
257
    this->set_temperature(temps[i], i);
9,669✔
258
  }
259

260
  // Read densities
261
  if (object_exists(cell_group, "density")) {
253✔
262
    vector<double> density;
198✔
263
    read_dataset(cell_group, "density", density);
198✔
264

265
    // Ensure number of densities makes sense
266
    auto n_density = density.size();
198!
267
    if (n_density > 1 && n_density != n_instances()) {
198!
268
      fatal_error(fmt::format("Number of densities for cell {} "
×
269
                              "doesn't match number of instances",
270
        id_));
×
271
    }
272

273
    // Set densities.
274
    for (int32_t i = 0; i < n_density; ++i) {
396✔
275
      this->set_density(density[i], i);
198✔
276
    }
277
  }
198✔
278

279
  close_group(cell_group);
253✔
280
}
253✔
281

282
void Cell::to_hdf5(hid_t cell_group) const
29,951✔
283
{
284

285
  // Create a group for this cell.
286
  auto group = create_group(cell_group, fmt::format("cell {}", id_));
29,951✔
287

288
  if (!name_.empty()) {
29,951✔
289
    write_string(group, "name", name_, false);
7,100✔
290
  }
291

292
  write_dataset(group, "universe", model::universes[universe_]->id_);
29,951✔
293

294
  to_hdf5_inner(group);
29,951✔
295

296
  // Write fill information.
297
  if (type_ == Fill::MATERIAL) {
29,951✔
298
    write_dataset(group, "fill_type", "material");
24,428✔
299
    std::vector<int32_t> mat_ids;
24,428✔
300
    for (auto i_mat : material_) {
50,149✔
301
      if (i_mat != MATERIAL_VOID) {
25,721✔
302
        mat_ids.push_back(model::materials[i_mat]->id_);
17,201✔
303
      } else {
304
        mat_ids.push_back(MATERIAL_VOID);
8,520✔
305
      }
306
    }
307
    if (mat_ids.size() == 1) {
24,428✔
308
      write_dataset(group, "material", mat_ids[0]);
24,237✔
309
    } else {
310
      write_dataset(group, "material", mat_ids);
191✔
311
    }
312

313
    std::vector<double> temps;
24,428✔
314
    for (auto sqrtkT_val : sqrtkT_)
60,502✔
315
      temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN);
36,074✔
316
    write_dataset(group, "temperature", temps);
24,428✔
317

318
    write_dataset(group, "density_mult", density_mult_);
24,428✔
319

320
  } else if (type_ == Fill::UNIVERSE) {
29,951✔
321
    write_dataset(group, "fill_type", "universe");
4,017✔
322
    write_dataset(group, "fill", model::universes[fill_]->id_);
4,017✔
323
    if (translation_ != Position(0, 0, 0)) {
4,017✔
324
      write_dataset(group, "translation", translation_);
1,837✔
325
    }
326
    if (!rotation_.empty()) {
4,017✔
327
      if (rotation_.size() == 12) {
264!
328
        std::array<double, 3> rot {rotation_[9], rotation_[10], rotation_[11]};
264✔
329
        write_dataset(group, "rotation", rot);
264✔
330
      } else {
331
        write_dataset(group, "rotation", rotation_);
×
332
      }
333
    }
334

335
  } else if (type_ == Fill::LATTICE) {
1,506!
336
    write_dataset(group, "fill_type", "lattice");
1,506✔
337
    write_dataset(group, "lattice", model::lattices[fill_]->id_);
1,506✔
338
  }
339

340
  close_group(group);
29,951✔
341
}
29,951✔
342

343
//==============================================================================
344
// XML parsing helpers for <cell> nodes
345
//==============================================================================
346

347
vector<int32_t> parse_cell_material_xml(pugi::xml_node node, int32_t cell_id)
28,468✔
348
{
349
  vector<std::string> mats {
28,468✔
350
    get_node_array<std::string>(node, "material", true)};
28,468✔
351
  if (mats.empty()) {
28,468!
352
    fatal_error(fmt::format(
×
353
      "An empty material element was specified for cell {}", cell_id));
354
  }
355
  vector<int32_t> material;
28,468✔
356
  material.reserve(mats.size());
28,468✔
357
  for (const auto& mat : mats) {
58,265✔
358
    if (mat == "void") {
29,797✔
359
      material.push_back(MATERIAL_VOID);
8,895✔
360
    } else {
361
      material.push_back(std::stoi(mat));
20,902✔
362
    }
363
  }
364
  return material;
28,468✔
365
}
28,468✔
366

367
vector<double> parse_cell_temperature_xml(pugi::xml_node node, int32_t cell_id)
431✔
368
{
369
  auto temperatures = get_node_array<double>(node, "temperature");
431✔
370
  if (temperatures.empty()) {
431!
371
    fatal_error(fmt::format(
×
372
      "An empty temperature element was specified for cell {}", cell_id));
373
  }
374
  for (auto T : temperatures) {
1,942✔
375
    if (T < 0) {
1,511!
376
      fatal_error(fmt::format(
×
377
        "Cell {} was specified with a negative temperature", cell_id));
378
    }
379
  }
380
  return temperatures;
431✔
381
}
×
382

383
vector<double> parse_cell_density_xml(pugi::xml_node node, int32_t cell_id)
75✔
384
{
385
  auto densities = get_node_array<double>(node, "density");
75✔
386
  if (densities.empty()) {
75!
387
    fatal_error(fmt::format(
×
388
      "An empty density element was specified for cell {}", cell_id));
389
  }
390
  for (auto rho : densities) {
1,230✔
391
    if (rho <= 0) {
1,155!
392
      fatal_error(fmt::format(
×
393
        "Cell {} was specified with a density less than or equal to zero",
394
        cell_id));
395
    }
396
  }
397
  return densities;
75✔
398
}
×
399

400
//==============================================================================
401
// CSGCell implementation
402
//==============================================================================
403

404
CSGCell::CSGCell(pugi::xml_node cell_node)
35,681✔
405
{
406
  if (check_for_node(cell_node, "id")) {
35,681!
407
    id_ = std::stoi(get_node_value(cell_node, "id"));
71,362✔
408
  } else {
409
    fatal_error("Must specify id of cell in geometry XML file.");
×
410
  }
411

412
  if (check_for_node(cell_node, "name")) {
35,681✔
413
    name_ = get_node_value(cell_node, "name");
9,075✔
414
  }
415

416
  if (check_for_node(cell_node, "universe")) {
35,681✔
417
    universe_ = std::stoi(get_node_value(cell_node, "universe"));
69,032✔
418
  } else {
419
    universe_ = 0;
1,165✔
420
  }
421

422
  // Make sure that either material or fill was specified, but not both.
423
  bool fill_present = check_for_node(cell_node, "fill");
35,681✔
424
  bool material_present = check_for_node(cell_node, "material");
35,681✔
425
  if (!(fill_present || material_present)) {
35,681!
426
    fatal_error(
×
427
      fmt::format("Neither material nor fill was specified for cell {}", id_));
×
428
  }
429
  if (fill_present && material_present) {
35,681!
430
    fatal_error(fmt::format("Cell {} has both a material and a fill specified; "
×
431
                            "only one can be specified per cell",
432
      id_));
×
433
  }
434

435
  if (fill_present) {
35,681✔
436
    fill_ = std::stoi(get_node_value(cell_node, "fill"));
14,444✔
437
    if (fill_ == universe_) {
7,222!
438
      fatal_error(fmt::format("Cell {} is filled with the same universe that "
×
439
                              "it is contained in.",
440
        id_));
×
441
    }
442
  } else {
443
    fill_ = C_NONE;
28,459✔
444
  }
445

446
  // Read the material element.  There can be zero materials (filled with a
447
  // universe), more than one material (distribmats), and some materials may
448
  // be "void".
449
  if (material_present) {
35,681✔
450
    material_ = parse_cell_material_xml(cell_node, id_);
28,459✔
451
  }
452

453
  // Read the temperature element which may be distributed like materials.
454
  if (check_for_node(cell_node, "temperature")) {
35,681✔
455
    sqrtkT_ = parse_cell_temperature_xml(cell_node, id_);
431✔
456
    sqrtkT_.shrink_to_fit();
431✔
457

458
    // Make sure this is a material-filled cell.
459
    if (material_.size() == 0) {
431!
460
      fatal_error(fmt::format(
×
461
        "Cell {} was specified with a temperature but no material. Temperature"
462
        "specification is only valid for cells filled with a material.",
463
        id_));
×
464
    }
465

466
    // Convert to sqrt(k*T).
467
    for (auto& T : sqrtkT_) {
1,942✔
468
      T = std::sqrt(K_BOLTZMANN * T);
1,511✔
469
    }
470
  }
471

472
  // Read the density element which can be distributed similar to temperature.
473
  // These get assigned to the density multiplier, requiring a division by
474
  // the material density.
475
  // Note: calculating the actual density multiplier is deferred until materials
476
  // are finalized. density_mult_ contains the true density in the meantime.
477
  if (check_for_node(cell_node, "density")) {
35,681✔
478
    density_mult_ = parse_cell_density_xml(cell_node, id_);
75✔
479
    density_mult_.shrink_to_fit();
75✔
480

481
    // Make sure this is a material-filled cell.
482
    if (material_.size() == 0) {
75!
483
      fatal_error(fmt::format(
×
484
        "Cell {} was specified with a density but no material. Density"
485
        "specification is only valid for cells filled with a material.",
486
        id_));
×
487
    }
488

489
    // Make sure this is a non-void material.
490
    for (auto mat_id : material_) {
150✔
491
      if (mat_id == MATERIAL_VOID) {
75!
492
        fatal_error(fmt::format(
×
493
          "Cell {} was specified with a density, but contains a void "
494
          "material. Density specification is only valid for cells "
495
          "filled with a non-void material.",
496
          id_));
×
497
      }
498
    }
499
  }
500

501
  // Read the region specification.
502
  std::string region_spec;
35,681✔
503
  if (check_for_node(cell_node, "region")) {
35,681✔
504
    region_spec = get_node_value(cell_node, "region");
26,582✔
505
  }
506

507
  // Get a tokenized representation of the region specification and apply De
508
  // Morgans law
509
  Region region(region_spec, id_);
35,681✔
510
  region_ = region;
35,681✔
511

512
  // Read the translation vector.
513
  if (check_for_node(cell_node, "translation")) {
35,681✔
514
    if (fill_ == C_NONE) {
2,557!
515
      fatal_error(fmt::format("Cannot apply a translation to cell {}"
×
516
                              " because it is not filled with another universe",
517
        id_));
×
518
    }
519

520
    auto xyz {get_node_array<double>(cell_node, "translation")};
2,557✔
521
    if (xyz.size() != 3) {
2,557!
522
      fatal_error(
×
523
        fmt::format("Non-3D translation vector applied to cell {}", id_));
×
524
    }
525
    translation_ = xyz;
2,557✔
526
  }
2,557✔
527

528
  // Read the rotation transform.
529
  if (check_for_node(cell_node, "rotation")) {
35,681✔
530
    auto rot {get_node_array<double>(cell_node, "rotation")};
390✔
531
    set_rotation(rot);
390✔
532
  }
390✔
533
}
35,681✔
534

535
//==============================================================================
536

537
void CSGCell::to_hdf5_inner(hid_t group_id) const
29,802✔
538
{
539
  write_string(group_id, "geom_type", "csg", false);
29,802✔
540
  write_string(group_id, "region", region_.str(), false);
29,802✔
541
}
29,802✔
542

543
//==============================================================================
544

545
vector<int32_t>::iterator CSGCell::find_left_parenthesis(
×
546
  vector<int32_t>::iterator start, const vector<int32_t>& infix)
547
{
548
  // start search at zero
549
  int parenthesis_level = 0;
×
550
  auto it = start;
×
551
  while (it != infix.begin()) {
×
552
    // look at two tokens at a time
553
    int32_t one = *it;
×
554
    int32_t two = *(it - 1);
×
555

556
    // decrement parenthesis level if there are two adjacent surfaces
557
    if (one < OP_UNION && two < OP_UNION) {
×
558
      parenthesis_level--;
×
559
      // increment if there are two adjacent operators
560
    } else if (one >= OP_UNION && two >= OP_UNION) {
×
561
      parenthesis_level++;
×
562
    }
563

564
    // if the level gets to zero, return the position
565
    if (parenthesis_level == 0) {
×
566
      // move the iterator back one before leaving the loop
567
      // so that all tokens in the parenthesis block are included
568
      it--;
×
569
      break;
570
    }
571

572
    // continue loop, one token at a time
573
    it--;
574
  }
575
  return it;
×
576
}
577

578
//==============================================================================
579
// Region implementation
580
//==============================================================================
581

582
Region::Region(std::string region_spec, int32_t cell_id)
35,780✔
583
{
584
  // Check if region_spec is not empty.
585
  if (!region_spec.empty()) {
35,780✔
586
    // Parse all halfspaces and operators except for intersection (whitespace).
587
    for (int i = 0; i < region_spec.size();) {
158,868✔
588
      if (region_spec[i] == '(') {
132,187✔
589
        expression_.push_back(OP_LEFT_PAREN);
1,681✔
590
        i++;
1,681✔
591

592
      } else if (region_spec[i] == ')') {
130,506✔
593
        expression_.push_back(OP_RIGHT_PAREN);
1,681✔
594
        i++;
1,681✔
595

596
      } else if (region_spec[i] == '|') {
128,825✔
597
        expression_.push_back(OP_UNION);
3,988✔
598
        i++;
3,988✔
599

600
      } else if (region_spec[i] == '~') {
124,837✔
601
        expression_.push_back(OP_COMPLEMENT);
30✔
602
        i++;
30✔
603

604
      } else if (region_spec[i] == '-' || region_spec[i] == '+' ||
210,682!
605
                 std::isdigit(region_spec[i])) {
85,875✔
606
        // This is the start of a halfspace specification.  Iterate j until we
607
        // find the end, then push-back everything between i and j.
608
        int j = i + 1;
73,004✔
609
        while (j < region_spec.size() && std::isdigit(region_spec[j])) {
148,661✔
610
          j++;
75,657✔
611
        }
612
        expression_.push_back(std::stoi(region_spec.substr(i, j - i)));
146,008✔
613
        i = j;
73,004✔
614

615
      } else if (std::isspace(region_spec[i])) {
51,803!
616
        i++;
51,803✔
617

618
      } else {
619
        auto err_msg =
×
620
          fmt::format("Region specification contains invalid character, \"{}\"",
621
            region_spec[i]);
×
622
        fatal_error(err_msg);
×
623
      }
×
624
    }
625

626
    // Add in intersection operators where a missing operator is needed.
627
    int i = 0;
628
    while (i < expression_.size() - 1) {
122,719✔
629
      bool left_compat {
96,038✔
630
        (expression_[i] < OP_UNION) || (expression_[i] == OP_RIGHT_PAREN)};
96,038✔
631
      bool right_compat {(expression_[i + 1] < OP_UNION) ||
96,038✔
632
                         (expression_[i + 1] == OP_LEFT_PAREN) ||
96,038✔
633
                         (expression_[i + 1] == OP_COMPLEMENT)};
5,729✔
634
      if (left_compat && right_compat) {
96,038✔
635
        expression_.insert(expression_.begin() + i + 1, OP_INTERSECTION);
42,335✔
636
      }
637
      i++;
638
    }
639

640
    // Remove complement operators using DeMorgan's laws
641
    auto it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT);
26,681✔
642
    while (it != expression_.end()) {
26,711✔
643
      // Erase complement
644
      expression_.erase(it);
30✔
645

646
      // Define stop given left parenthesis or not
647
      auto stop = it;
30✔
648
      if (*it == OP_LEFT_PAREN) {
30!
649
        int depth = 1;
650
        do {
240✔
651
          stop++;
240✔
652
          if (*stop > OP_COMPLEMENT) {
240✔
653
            if (*stop == OP_RIGHT_PAREN) {
30!
654
              depth--;
30✔
655
            } else {
656
              depth++;
×
657
            }
658
          }
659
        } while (depth > 0);
240✔
660
        it++;
30✔
661
      }
662

663
      // apply DeMorgan's law to any surfaces/operators between these
664
      // positions in the RPN
665
      apply_demorgan(it, stop);
30✔
666
      // update iterator position
667
      it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT);
30✔
668
    }
669

670
    // Convert user IDs to surface indices.
671
    for (auto& r : expression_) {
149,370✔
672
      if (r < OP_UNION) {
122,689✔
673
        const auto& it {model::surface_map.find(abs(r))};
73,004!
674
        if (it == model::surface_map.end()) {
73,004!
675
          throw std::runtime_error {
×
676
            "Invalid surface ID " + std::to_string(abs(r)) +
×
677
            " specified in region for cell " + std::to_string(cell_id) + "."};
×
678
        }
679
        r = (r > 0) ? it->second + 1 : -(it->second + 1);
73,004✔
680
      }
681
    }
682

683
    // Check if this is a simple cell.
684
    simple_ = true;
26,681✔
685
    for (int32_t token : expression_) {
134,149✔
686
      if (token == OP_UNION) {
108,677✔
687
        simple_ = false;
1,209✔
688
        // Ensure intersections have precedence over unions
689
        enforce_precedence();
1,209✔
690
        break;
691
      }
692
    }
693

694
    // If this cell is simple, remove all the superfluous operator tokens.
695
    if (simple_) {
26,681✔
696
      for (auto it = expression_.begin(); it != expression_.end(); it++) {
126,238✔
697
        if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) {
100,766!
698
          expression_.erase(it);
37,647✔
699
          it--;
100,766✔
700
        }
701
      }
702
    }
703
    expression_.shrink_to_fit();
26,681✔
704

705
  } else {
706
    simple_ = true;
9,099✔
707
  }
708
}
35,780✔
709

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

712
void Region::apply_demorgan(
30✔
713
  vector<int32_t>::iterator start, vector<int32_t>::iterator stop)
714
{
715
  do {
210✔
716
    if (*start < OP_UNION) {
210✔
717
      *start *= -1;
120✔
718
    } else if (*start == OP_UNION) {
90!
719
      *start = OP_INTERSECTION;
×
720
    } else if (*start == OP_INTERSECTION) {
90!
721
      *start = OP_UNION;
90✔
722
    }
723
    start++;
210✔
724
  } while (start < stop);
210✔
725
}
30✔
726

727
//==============================================================================
728
//! Add precedence for infix regions so intersections have higher
729
//! precedence than unions using parentheses.
730
//==============================================================================
731

732
void Region::add_parentheses(int64_t start)
96✔
733
{
734
  int32_t start_token = expression_[start];
96!
735
  // Add left parenthesis and set new position to be after parenthesis
736
  if (start_token == OP_UNION) {
96!
737
    start += 2;
×
738
  }
739
  expression_.insert(expression_.begin() + start - 1, OP_LEFT_PAREN);
96✔
740

741
  // Add right parenthesis
742
  // While the start iterator is within the bounds of infix
743
  while (start + 1 < expression_.size()) {
430✔
744
    start++;
408✔
745

746
    // If the current token is an operator and is different than the start token
747
    if (expression_[start] >= OP_UNION && expression_[start] != start_token) {
408✔
748
      // Skip wrapped regions but save iterator position to check precedence and
749
      // add right parenthesis, right parenthesis position depends on the
750
      // operator, when the operator is a union then do not include the operator
751
      // in the region, when the operator is an intersection then include the
752
      // operator and next surface
753
      if (expression_[start] == OP_LEFT_PAREN) {
85✔
754
        int depth = 1;
755
        do {
44✔
756
          start++;
44✔
757
          if (expression_[start] > OP_COMPLEMENT) {
44✔
758
            if (expression_[start] == OP_RIGHT_PAREN) {
11!
759
              depth--;
11✔
760
            } else {
761
              depth++;
×
762
            }
763
          }
764
        } while (depth > 0);
44✔
765
      } else {
766
        if (start_token == OP_UNION) {
74!
767
          --start;
×
768
        }
769
        expression_.insert(expression_.begin() + start, OP_RIGHT_PAREN);
74✔
770
        return;
74✔
771
      }
772
    }
773
  }
774
  // If we get here a right parenthesis hasn't been placed
775
  expression_.push_back(OP_RIGHT_PAREN);
22✔
776
}
777

778
//==============================================================================
779
//! Add parentheses to enforce operator precedence in region expressions
780
//!
781
//! This function ensures that intersection operators have higher precedence
782
//! than union operators by adding parentheses where needed. For example:
783
//!   "1 2 | 3" becomes "(1 2) | 3"
784
//!   "1 | 2 3" becomes "1 | (2 3)"
785
//!
786
//! The algorithm uses stacks to track the current operator type and its
787
//! position at each parenthesis depth level. When it encounters a different
788
//! operator at the same depth, it adds parentheses to group the
789
//! higher-precedence operations.
790
//==============================================================================
791

792
void Region::enforce_precedence()
1,209✔
793
{
794
  // Stack tracking the operator type at each depth (0 = no operator seen yet)
795
  vector<int32_t> op_stack = {0};
1,209✔
796

797
  // Stack tracking where the operator sequence started at each depth
798
  vector<std::size_t> pos_stack = {0};
1,209✔
799

800
  for (int64_t i = 0; i < expression_.size(); ++i) {
24,149✔
801
    int32_t token = expression_[i];
22,940✔
802

803
    if (token == OP_LEFT_PAREN) {
22,940✔
804
      // Entering a new parenthesis level - push new tracking state
805
      op_stack.push_back(0);
1,877✔
806
      pos_stack.push_back(0);
1,877✔
807
      continue;
1,877✔
808
    } else if (token == OP_RIGHT_PAREN) {
21,063✔
809
      // Exiting a parenthesis level - pop tracking state (keep at least one)
810
      if (op_stack.size() > 1) {
1,836!
811
        op_stack.pop_back();
1,836✔
812
        pos_stack.pop_back();
1,836✔
813
      }
814
      continue;
1,836✔
815
    }
816

817
    if (token == OP_UNION || token == OP_INTERSECTION) {
19,227✔
818
      if (op_stack.back() == 0) {
9,009✔
819
        // First operator at this depth - record it and its position
820
        op_stack.back() = token;
3,130✔
821
        pos_stack.back() = i;
3,130✔
822
      } else if (token != op_stack.back()) {
5,879✔
823
        // Encountered a different operator at the same depth - need to add
824
        // parentheses to enforce precedence. Intersection has higher
825
        // precedence, so we parenthesize the intersection terms.
826
        if (op_stack.back() == OP_INTERSECTION) {
96✔
827
          add_parentheses(pos_stack.back());
48✔
828
        } else {
829
          add_parentheses(i);
48✔
830
        }
831

832
        // Restart the scan since we modified the expression
833
        i = -1; // Will be incremented to 0 by the for loop
96✔
834
        op_stack = {0};
96✔
835
        pos_stack = {0};
96✔
836
      }
837
    }
838
  }
839
}
1,209✔
840

841
//==============================================================================
842
//! Convert infix region specification to Reverse Polish Notation (RPN)
843
//!
844
//! This function uses the shunting-yard algorithm.
845
//==============================================================================
846

847
vector<int32_t> Region::generate_postfix(int32_t cell_id) const
44✔
848
{
849
  vector<int32_t> rpn;
44✔
850
  vector<int32_t> stack;
44✔
851

852
  for (int32_t token : expression_) {
990✔
853
    if (token < OP_UNION) {
946✔
854
      // If token is not an operator, add it to output
855
      rpn.push_back(token);
396✔
856
    } else if (token < OP_RIGHT_PAREN) {
550✔
857
      // Regular operators union, intersection, complement
858
      while (stack.size() > 0) {
561✔
859
        int32_t op = stack.back();
462✔
860

861
        if (op < OP_RIGHT_PAREN && ((token == OP_COMPLEMENT && token < op) ||
462!
862
                                     (token != OP_COMPLEMENT && token <= op))) {
209!
863
          // While there is an operator, op, on top of the stack, if the token
864
          // is left-associative and its precedence is less than or equal to
865
          // that of op or if the token is right-associative and its precedence
866
          // is less than that of op, move op to the output queue and push the
867
          // token on to the stack. Note that only complement is
868
          // right-associative.
869
          rpn.push_back(op);
209✔
870
          stack.pop_back();
209✔
871
        } else {
872
          break;
873
        }
874
      }
875

876
      stack.push_back(token);
352✔
877

878
    } else if (token == OP_LEFT_PAREN) {
198✔
879
      // If the token is a left parenthesis, push it onto the stack
880
      stack.push_back(token);
99✔
881

882
    } else {
883
      // If the token is a right parenthesis, move operators from the stack to
884
      // the output queue until reaching the left parenthesis.
885
      for (auto it = stack.rbegin(); *it != OP_LEFT_PAREN; it++) {
198✔
886
        // If we run out of operators without finding a left parenthesis, it
887
        // means there are mismatched parentheses.
888
        if (it == stack.rend()) {
99!
889
          fatal_error(fmt::format(
×
890
            "Mismatched parentheses in region specification for cell {}",
891
            cell_id));
892
        }
893
        rpn.push_back(stack.back());
99✔
894
        stack.pop_back();
99✔
895
      }
896

897
      // Pop the left parenthesis.
898
      stack.pop_back();
946✔
899
    }
900
  }
901

902
  while (stack.size() > 0) {
44✔
903
    int32_t op = stack.back();
44!
904

905
    // If the operator is a parenthesis it is mismatched.
906
    if (op >= OP_RIGHT_PAREN) {
44!
907
      fatal_error(fmt::format(
×
908
        "Mismatched parentheses in region specification for cell {}", cell_id));
909
    }
910

911
    rpn.push_back(stack.back());
44✔
912
    stack.pop_back();
88✔
913
  }
914

915
  return rpn;
44✔
916
}
44✔
917

918
//==============================================================================
919

920
std::string Region::str() const
29,901✔
921
{
922
  std::stringstream region_spec {};
29,901✔
923
  if (!expression_.empty()) {
29,901✔
924
    for (int32_t token : expression_) {
91,501✔
925
      if (token == OP_LEFT_PAREN) {
69,984✔
926
        region_spec << " (";
1,505✔
927
      } else if (token == OP_RIGHT_PAREN) {
68,479✔
928
        region_spec << " )";
1,505✔
929
      } else if (token == OP_COMPLEMENT) {
66,974!
930
        region_spec << " ~";
×
931
      } else if (token == OP_INTERSECTION) {
66,974✔
932
      } else if (token == OP_UNION) {
63,132✔
933
        region_spec << " |";
3,534✔
934
      } else {
935
        // Note the off-by-one indexing
936
        auto surf_id = model::surfaces[abs(token) - 1]->id_;
59,598✔
937
        region_spec << " " << ((token > 0) ? surf_id : -surf_id);
59,598✔
938
      }
939
    }
940
  }
941
  return region_spec.str();
59,802✔
942
}
29,901✔
943

944
//==============================================================================
945

946
std::pair<double, int32_t> Region::distance(
2,147,483,647✔
947
  Position r, Direction u, int32_t on_surface) const
948
{
949
  if (simple_) {
2,147,483,647✔
950
    return distance_simple(r, u, on_surface);
2,147,483,647✔
951
  } else {
952
    return distance_complex(r, u, on_surface);
18,872,603✔
953
  }
954
}
955

956
//==============================================================================
957

958
std::pair<double, int32_t> Region::distance_simple(
2,147,483,647✔
959
  Position r, Direction u, int32_t on_surface) const
960
{
961
  double min_dist {INFTY};
2,147,483,647✔
962
  int32_t i_surf {std::numeric_limits<int32_t>::max()};
2,147,483,647✔
963

964
  for (int32_t token : expression_) {
2,147,483,647✔
965
    // Ignore this token if it corresponds to an operator rather than a region.
966
    if (token >= OP_UNION)
2,147,483,647!
UNCOV
967
      continue;
×
968

969
    // Calculate the distance to this surface.
970
    // Note the off-by-one indexing
971
    bool coincident {std::abs(token) == std::abs(on_surface)};
2,147,483,647✔
972
    double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
2,147,483,647✔
973

974
    // Check if this distance is the new minimum.
975
    if (d < min_dist) {
2,147,483,647✔
976
      if (min_dist - d >= FP_PRECISION * min_dist) {
2,147,483,647!
977
        min_dist = d;
2,147,483,647✔
978
        i_surf = -token;
2,147,483,647✔
979
      }
980
    }
981
  }
982

983
  return {min_dist, i_surf};
2,147,483,647✔
984
}
985

986
//==============================================================================
987

988
std::pair<double, int32_t> Region::distance_complex(
18,872,603✔
989
  Position r, Direction u, int32_t on_surface) const
990
{
991
  double min_dist;
18,872,603✔
992
  int32_t i_surf;
18,872,603✔
993
  double atleast {-1.0};
18,872,603✔
994
  bool in_region = contains_complex(r, u, on_surface);
18,872,603✔
995

996
  while (true) {
46,385,923✔
997
    min_dist = INFTY;
32,629,263✔
998
    i_surf = std::numeric_limits<int32_t>::max();
32,629,263✔
999

1000
    for (int32_t token : expression_) {
1,358,209,338✔
1001
      // Ignore this token if it corresponds to an operator rather than a
1002
      // region.
1003
      if (token >= OP_UNION)
1,325,580,075✔
1004
        continue;
794,831,939✔
1005

1006
      // Calculate the distance to this surface.
1007
      // Note the off-by-one indexing
1008
      bool coincident {std::abs(token) == std::abs(on_surface)};
530,748,136✔
1009
      double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
530,748,136✔
1010

1011
      // Check if this distance is the new minimum.
1012
      if ((d > atleast) && (d < min_dist)) {
530,748,136✔
1013
        if (min_dist - d >= FP_PRECISION * min_dist) {
66,506,896!
1014
          min_dist = d;
66,506,896✔
1015
          i_surf = -token;
66,506,896✔
1016
        }
1017
      }
1018
    }
1019
    if (min_dist == INFTY)
32,629,263✔
1020
      break;
1021
    auto p = r + (min_dist + TINY_BIT) * u;
30,884,045✔
1022
    if (contains_complex(p, u, on_surface) != in_region)
30,884,045✔
1023
      break;
1024
    atleast = min_dist;
13,756,660✔
1025
  }
13,756,660✔
1026
  return {min_dist, i_surf};
18,872,603✔
1027
}
1028

1029
//==============================================================================
1030

1031
bool Region::contains(Position r, Direction u, int32_t on_surface) const
2,147,483,647✔
1032
{
1033
  if (simple_) {
2,147,483,647✔
1034
    return contains_simple(r, u, on_surface);
2,147,483,647✔
1035
  } else {
1036
    return contains_complex(r, u, on_surface);
14,026,985✔
1037
  }
1038
}
1039

1040
//==============================================================================
1041

1042
bool Region::contains_simple(Position r, Direction u, int32_t on_surface) const
2,147,483,647✔
1043
{
1044
  for (int32_t token : expression_) {
2,147,483,647✔
1045
    // Assume that no tokens are operators. Evaluate the sense of particle with
1046
    // respect to the surface and see if the token matches the sense. If the
1047
    // particle's surface attribute is set and matches the token, that
1048
    // overrides the determination based on sense().
1049
    if (token == on_surface) {
2,147,483,647✔
1050
    } else if (-token == on_surface) {
2,147,483,647✔
1051
      return false;
1052
    } else {
1053
      // Note the off-by-one indexing
1054
      bool sense = model::surfaces[abs(token) - 1]->sense(r, u);
2,147,483,647✔
1055
      if (sense != (token > 0)) {
2,147,483,647✔
1056
        return false;
1057
      }
1058
    }
1059
  }
1060
  return true;
1061
}
1062

1063
//==============================================================================
1064

1065
bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const
63,783,633✔
1066
{
1067
  bool in_cell = true;
63,783,633✔
1068
  int total_depth = 0;
63,783,633✔
1069

1070
  // For each token
1071
  for (auto it = expression_.begin(); it != expression_.end(); it++) {
1,141,248,788✔
1072
    int32_t token = *it;
1,099,159,324✔
1073

1074
    // If the token is a surface evaluate the sense
1075
    // If the token is a union or intersection check to
1076
    // short circuit
1077
    if (token < OP_UNION) {
1,099,159,324✔
1078
      if (token == on_surface) {
425,418,552✔
1079
        in_cell = true;
1080
      } else if (-token == on_surface) {
398,701,904✔
1081
        in_cell = false;
1082
      } else {
1083
        // Note the off-by-one indexing
1084
        bool sense = model::surfaces[abs(token) - 1]->sense(r, u);
395,672,051✔
1085
        in_cell = (sense == (token > 0));
395,672,051✔
1086
      }
1087
    } else if ((token == OP_UNION && in_cell == true) ||
673,740,772✔
1088
               (token == OP_INTERSECTION && in_cell == false)) {
326,760,122✔
1089
      // If the total depth is zero return
1090
      if (total_depth == 0) {
110,408,595✔
1091
        return in_cell;
21,694,169✔
1092
      }
1093

1094
      total_depth--;
88,714,426✔
1095

1096
      // While the iterator is within the bounds of the vector
1097
      int depth = 1;
88,714,426✔
1098
      do {
1,375,924,514✔
1099
        // Get next token
1100
        it++;
1,375,924,514✔
1101
        int32_t next_token = *it;
1,375,924,514✔
1102

1103
        // If the token is an a parenthesis
1104
        if (next_token > OP_COMPLEMENT) {
1,375,924,514✔
1105
          // Adjust depth accordingly
1106
          if (next_token == OP_RIGHT_PAREN) {
390,386,688✔
1107
            depth--;
239,550,557✔
1108
          } else {
1109
            depth++;
150,836,131✔
1110
          }
1111
        }
1112
      } while (depth > 0);
1,375,924,514✔
1113
    } else if (token == OP_LEFT_PAREN) {
563,332,177✔
1114
      total_depth++;
145,205,842✔
1115
    } else if (token == OP_RIGHT_PAREN) {
418,126,335✔
1116
      total_depth--;
56,491,416✔
1117
    }
1118
  }
1119
  return in_cell;
1120
}
1121

1122
//==============================================================================
1123

1124
BoundingBox Region::bounding_box(int32_t cell_id) const
88✔
1125
{
1126
  if (simple_) {
88✔
1127
    return bounding_box_simple();
44✔
1128
  } else {
1129
    auto postfix = generate_postfix(cell_id);
44✔
1130
    return bounding_box_complex(postfix);
88✔
1131
  }
44✔
1132
}
1133

1134
//==============================================================================
1135

1136
BoundingBox Region::bounding_box_simple() const
44✔
1137
{
1138
  BoundingBox bbox;
44✔
1139
  for (int32_t token : expression_) {
176✔
1140
    bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0);
132✔
1141
  }
1142
  return bbox;
44✔
1143
}
1144

1145
//==============================================================================
1146

1147
BoundingBox Region::bounding_box_complex(vector<int32_t> postfix) const
44✔
1148
{
1149
  vector<BoundingBox> stack(postfix.size());
44✔
1150
  int i_stack = -1;
44✔
1151

1152
  for (auto& token : postfix) {
792✔
1153
    if (token == OP_UNION) {
748✔
1154
      stack[i_stack - 1] = stack[i_stack - 1] | stack[i_stack];
154✔
1155
      i_stack--;
154✔
1156
    } else if (token == OP_INTERSECTION) {
594✔
1157
      stack[i_stack - 1] = stack[i_stack - 1] & stack[i_stack];
198✔
1158
      i_stack--;
198✔
1159
    } else {
1160
      i_stack++;
396✔
1161
      stack[i_stack] = model::surfaces[abs(token) - 1]->bounding_box(token > 0);
396✔
1162
    }
1163
  }
1164

1165
  assert(i_stack == 0);
44!
1166
  return stack.front();
44✔
1167
}
44✔
1168

1169
//==============================================================================
1170

1171
vector<int32_t> Region::surfaces() const
5,270✔
1172
{
1173
  if (simple_) {
5,270✔
1174
    return expression_;
5,250✔
1175
  }
1176

1177
  vector<int32_t> surfaces = expression_;
20✔
1178

1179
  auto it = std::find_if(surfaces.begin(), surfaces.end(),
20✔
1180
    [&](const auto& value) { return value >= OP_UNION; });
20!
1181

1182
  while (it != surfaces.end()) {
60✔
1183
    surfaces.erase(it);
40✔
1184

1185
    it = std::find_if(surfaces.begin(), surfaces.end(),
40✔
1186
      [&](const auto& value) { return value >= OP_UNION; });
80!
1187
  }
1188

1189
  return surfaces;
20✔
1190
}
5,270✔
1191

1192
//==============================================================================
1193
// Non-method functions
1194
//==============================================================================
1195

1196
void read_cells(pugi::xml_node node)
8,739✔
1197
{
1198
  // Count the number of cells.
1199
  int n_cells = 0;
8,739✔
1200
  for (pugi::xml_node cell_node : node.children("cell")) {
44,420✔
1201
    n_cells++;
35,681✔
1202
  }
1203

1204
  // Loop over XML cell elements and populate the array.
1205
  model::cells.reserve(n_cells);
8,739✔
1206
  for (pugi::xml_node cell_node : node.children("cell")) {
44,420✔
1207
    model::cells.push_back(make_unique<CSGCell>(cell_node));
35,681✔
1208
  }
1209

1210
  // Fill the cell map.
1211
  for (int i = 0; i < model::cells.size(); i++) {
44,420✔
1212
    int32_t id = model::cells[i]->id_;
35,681!
1213
    auto search = model::cell_map.find(id);
35,681!
1214
    if (search == model::cell_map.end()) {
35,681!
1215
      model::cell_map[id] = i;
35,681✔
1216
    } else {
1217
      fatal_error(
×
1218
        fmt::format("Two or more cells use the same unique ID: {}", id));
×
1219
    }
1220
  }
1221

1222
  read_dagmc_universes(node);
8,739✔
1223

1224
  populate_universes();
8,737✔
1225

1226
  // Allocate the cell overlap count if necessary.
1227
  if (settings::check_overlaps) {
8,737✔
1228
    model::overlap_check_count.resize(model::cells.size(), 0);
119✔
1229
  }
1230

1231
  if (model::cells.size() == 0) {
8,737!
1232
    fatal_error("No cells were found in the geometry.xml file");
×
1233
  }
1234
}
8,737✔
1235

1236
void populate_universes()
8,739✔
1237
{
1238
  // Used to map universe index to the index of an implicit complement cell for
1239
  // DAGMC universes
1240
  std::unordered_map<int, int> implicit_comp_cells;
8,739✔
1241

1242
  // Populate the Universe vector and map.
1243
  for (int index_cell = 0; index_cell < model::cells.size(); index_cell++) {
44,621✔
1244
    int32_t uid = model::cells[index_cell]->universe_;
35,882✔
1245
    auto it = model::universe_map.find(uid);
35,882✔
1246
    if (it == model::universe_map.end()) {
35,882✔
1247
      model::universes.push_back(make_unique<Universe>());
40,990✔
1248
      model::universes.back()->id_ = uid;
20,495✔
1249
      model::universes.back()->cells_.push_back(index_cell);
20,495✔
1250
      model::universe_map[uid] = model::universes.size() - 1;
20,495✔
1251
    } else {
1252
#ifdef OPENMC_DAGMC_ENABLED
1253
      // Skip implicit complement cells for now
1254
      Universe* univ = model::universes[it->second].get();
2,044!
1255
      DAGUniverse* dag_univ = dynamic_cast<DAGUniverse*>(univ);
2,044!
1256
      if (dag_univ && (dag_univ->implicit_complement_idx() == index_cell)) {
2,044✔
1257
        implicit_comp_cells[it->second] = index_cell;
43✔
1258
        continue;
43✔
1259
      }
1260
#endif
1261

1262
      model::universes[it->second]->cells_.push_back(index_cell);
15,344✔
1263
    }
1264
  }
1265

1266
  // Add DAGUniverse implicit complement cells last
1267
  for (const auto& it : implicit_comp_cells) {
8,782✔
1268
    int index_univ = it.first;
43✔
1269
    int index_cell = it.second;
43✔
1270
    model::universes[index_univ]->cells_.push_back(index_cell);
43!
1271
  }
1272

1273
  model::universes.shrink_to_fit();
8,739✔
1274
}
8,739✔
1275

1276
//==============================================================================
1277
// C-API functions
1278
//==============================================================================
1279

1280
extern "C" int openmc_cell_get_fill(
220✔
1281
  int32_t index, int* type, int32_t** indices, int32_t* n)
1282
{
1283
  if (index >= 0 && index < model::cells.size()) {
220!
1284
    Cell& c {*model::cells[index]};
220✔
1285
    *type = static_cast<int>(c.type_);
220✔
1286
    if (c.type_ == Fill::MATERIAL) {
220✔
1287
      *indices = c.material_.data();
209✔
1288
      *n = c.material_.size();
209✔
1289
    } else {
1290
      *indices = &c.fill_;
11✔
1291
      *n = 1;
11✔
1292
    }
1293
  } else {
1294
    set_errmsg("Index in cells array is out of bounds.");
×
1295
    return OPENMC_E_OUT_OF_BOUNDS;
×
1296
  }
1297
  return 0;
1298
}
1299

1300
extern "C" int openmc_cell_set_fill(
11✔
1301
  int32_t index, int type, int32_t n, const int32_t* indices)
1302
{
1303
  Fill filltype = static_cast<Fill>(type);
11✔
1304
  if (index >= 0 && index < model::cells.size()) {
11!
1305
    Cell& c {*model::cells[index]};
11!
1306
    if (filltype == Fill::MATERIAL) {
11!
1307
      c.type_ = Fill::MATERIAL;
11✔
1308
      c.material_.clear();
11!
1309
      for (int i = 0; i < n; i++) {
22✔
1310
        int i_mat = indices[i];
11✔
1311
        if (i_mat == MATERIAL_VOID) {
11!
1312
          c.material_.push_back(MATERIAL_VOID);
×
1313
        } else if (i_mat >= 0 && i_mat < model::materials.size()) {
11!
1314
          c.material_.push_back(i_mat);
11✔
1315
        } else {
1316
          set_errmsg("Index in materials array is out of bounds.");
×
1317
          return OPENMC_E_OUT_OF_BOUNDS;
×
1318
        }
1319
      }
1320
      c.material_.shrink_to_fit();
11✔
1321
    } else if (filltype == Fill::UNIVERSE) {
×
1322
      c.type_ = Fill::UNIVERSE;
×
1323
    } else {
1324
      c.type_ = Fill::LATTICE;
×
1325
    }
1326
  } else {
1327
    set_errmsg("Index in cells array is out of bounds.");
×
1328
    return OPENMC_E_OUT_OF_BOUNDS;
×
1329
  }
1330
  return 0;
1331
}
1332

1333
extern "C" int openmc_cell_set_temperature(
88✔
1334
  int32_t index, double T, const int32_t* instance, bool set_contained)
1335
{
1336
  if (index < 0 || index >= model::cells.size()) {
88!
1337
    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
×
1338
    return OPENMC_E_OUT_OF_BOUNDS;
×
1339
  }
1340

1341
  int32_t instance_index = instance ? *instance : -1;
88✔
1342
  try {
88✔
1343
    model::cells[index]->set_temperature(T, instance_index, set_contained);
88✔
1344
  } catch (const std::exception& e) {
×
1345
    set_errmsg(e.what());
×
1346
    return OPENMC_E_UNASSIGNED;
×
1347
  }
×
1348
  return 0;
1349
}
1350

1351
extern "C" int openmc_cell_set_density(
88✔
1352
  int32_t index, double density, const int32_t* instance, bool set_contained)
1353
{
1354
  if (index < 0 || index >= model::cells.size()) {
88!
1355
    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
×
1356
    return OPENMC_E_OUT_OF_BOUNDS;
×
1357
  }
1358

1359
  int32_t instance_index = instance ? *instance : -1;
88✔
1360
  try {
88✔
1361
    model::cells[index]->set_density(density, instance_index, set_contained);
88✔
1362
  } catch (const std::exception& e) {
×
1363
    set_errmsg(e.what());
×
1364
    return OPENMC_E_UNASSIGNED;
×
1365
  }
×
1366
  return 0;
1367
}
1368

1369
extern "C" int openmc_cell_get_temperature(
9,628✔
1370
  int32_t index, const int32_t* instance, double* T)
1371
{
1372
  if (index < 0 || index >= model::cells.size()) {
9,628!
1373
    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
×
1374
    return OPENMC_E_OUT_OF_BOUNDS;
×
1375
  }
1376

1377
  int32_t instance_index = instance ? *instance : -1;
9,628✔
1378
  try {
9,628✔
1379
    *T = model::cells[index]->temperature(instance_index);
9,628✔
1380
  } catch (const std::exception& e) {
×
1381
    set_errmsg(e.what());
×
1382
    return OPENMC_E_UNASSIGNED;
×
1383
  }
×
1384
  return 0;
9,628✔
1385
}
1386

1387
extern "C" int openmc_cell_get_density(
88✔
1388
  int32_t index, const int32_t* instance, double* density)
1389
{
1390
  if (index < 0 || index >= model::cells.size()) {
88!
1391
    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
×
1392
    return OPENMC_E_OUT_OF_BOUNDS;
×
1393
  }
1394

1395
  int32_t instance_index = instance ? *instance : -1;
88✔
1396
  try {
88✔
1397
    if (model::cells[index]->type_ != Fill::MATERIAL) {
88!
1398
      fatal_error(
×
1399
        fmt::format("Cell {}, instance {} is not filled with a material.",
×
1400
          model::cells[index]->id_, instance_index));
×
1401
    }
1402

1403
    int32_t mat_index = model::cells[index]->material(instance_index);
88!
1404
    if (mat_index == MATERIAL_VOID) {
88!
1405
      *density = 0.0;
×
1406
    } else {
1407
      *density = model::cells[index]->density_mult(instance_index) *
88✔
1408
                 model::materials[mat_index]->density_gpcc();
176!
1409
    }
1410
  } catch (const std::exception& e) {
×
1411
    set_errmsg(e.what());
×
1412
    return OPENMC_E_UNASSIGNED;
×
1413
  }
×
1414
  return 0;
1415
}
1416

1417
//! Get the bounding box of a cell
1418
extern "C" int openmc_cell_bounding_box(
55✔
1419
  const int32_t index, double* llc, double* urc)
1420
{
1421

1422
  BoundingBox bbox;
55✔
1423

1424
  const auto& c = model::cells[index];
55✔
1425
  bbox = c->bounding_box();
55✔
1426

1427
  // set lower left corner values
1428
  llc[0] = bbox.min.x;
55✔
1429
  llc[1] = bbox.min.y;
55✔
1430
  llc[2] = bbox.min.z;
55✔
1431

1432
  // set upper right corner values
1433
  urc[0] = bbox.max.x;
55✔
1434
  urc[1] = bbox.max.y;
55✔
1435
  urc[2] = bbox.max.z;
55✔
1436

1437
  return 0;
55✔
1438
}
1439

1440
//! Get the name of a cell
1441
extern "C" int openmc_cell_get_name(int32_t index, const char** name)
413✔
1442
{
1443
  if (index < 0 || index >= model::cells.size()) {
413!
1444
    set_errmsg("Index in cells array is out of bounds.");
×
1445
    return OPENMC_E_OUT_OF_BOUNDS;
×
1446
  }
1447

1448
  *name = model::cells[index]->name().data();
413✔
1449

1450
  return 0;
413✔
1451
}
1452

1453
//! Set the name of a cell
1454
extern "C" int openmc_cell_set_name(int32_t index, const char* name)
11✔
1455
{
1456
  if (index < 0 || index >= model::cells.size()) {
11!
1457
    set_errmsg("Index in cells array is out of bounds.");
×
1458
    return OPENMC_E_OUT_OF_BOUNDS;
×
1459
  }
1460

1461
  model::cells[index]->set_name(name);
22✔
1462

1463
  return 0;
11✔
1464
}
1465

1466
//==============================================================================
1467
//! Define a containing (parent) cell
1468
//==============================================================================
1469

1470
//! Used to locate a universe fill in the geometry
1471
struct ParentCell {
1472
  bool operator==(const ParentCell& other) const
135✔
1473
  {
1474
    return cell_index == other.cell_index &&
135!
1475
           lattice_index == other.lattice_index;
135!
1476
  }
1477

1478
  bool operator<(const ParentCell& other) const
1479
  {
1480
    return cell_index < other.cell_index ||
1481
           (cell_index == other.cell_index &&
1482
             lattice_index < other.lattice_index);
1483
  }
1484

1485
  int64_t cell_index;
1486
  int64_t lattice_index;
1487
};
1488

1489
//! Structure used to insert ParentCell into hashed STL data structures
1490
struct ParentCellHash {
1491
  std::size_t operator()(const ParentCell& p) const
661✔
1492
  {
1493
    return 4096 * p.cell_index + p.lattice_index;
661!
1494
  }
1495
};
1496

1497
//! Used to manage a traversal stack when locating parent cells of a cell
1498
//! instance in the model
1499
struct ParentCellStack {
136✔
1500

1501
  //! push method that adds to the parent_cells visited cells for this search
1502
  //! universe
1503
  void push(int32_t search_universe, const ParentCell& pc)
105✔
1504
  {
1505
    parent_cells_.push_back(pc);
105✔
1506
    // add parent cell to the set of cells we've visited for this search
1507
    // universe
1508
    visited_cells_[search_universe].insert(pc);
105✔
1509
  }
105✔
1510

1511
  //! removes the last parent_cell and clears the visited cells for the popped
1512
  //! cell's universe
1513
  void pop()
75✔
1514
  {
1515
    visited_cells_[this->current_univ()].clear();
75✔
1516
    parent_cells_.pop_back();
75✔
1517
  }
75✔
1518

1519
  //! checks whether or not the parent cell has been visited already for this
1520
  //! search universe
1521
  bool visited(int32_t search_universe, const ParentCell& parent_cell)
556✔
1522
  {
1523
    return visited_cells_[search_universe].count(parent_cell) != 0;
556✔
1524
  }
1525

1526
  //! return the next universe to search for a parent cell
1527
  int32_t current_univ() const
75✔
1528
  {
1529
    return model::cells[parent_cells_.back().cell_index]->universe_;
75✔
1530
  }
1531

1532
  //! indicates whether nor not parent cells are present on the stack
1533
  bool empty() const { return parent_cells_.empty(); }
75✔
1534

1535
  //! compute an instance for the provided distribcell index
1536
  int32_t compute_instance(int32_t distribcell_index) const
211✔
1537
  {
1538
    if (distribcell_index == C_NONE)
211✔
1539
      return 0;
1540

1541
    int32_t instance = 0;
120✔
1542
    for (const auto& parent_cell : this->parent_cells_) {
225✔
1543
      auto& cell = model::cells[parent_cell.cell_index];
105!
1544
      if (cell->type_ == Fill::UNIVERSE) {
105!
1545
        instance += cell->offset_[distribcell_index];
×
1546
      } else if (cell->type_ == Fill::LATTICE) {
105!
1547
        auto& lattice = model::lattices[cell->fill_];
105✔
1548
        instance +=
105✔
1549
          lattice->offset(distribcell_index, parent_cell.lattice_index);
105✔
1550
      }
1551
    }
1552
    return instance;
1553
  }
1554

1555
  // Accessors
1556
  vector<ParentCell>& parent_cells() { return parent_cells_; }
136✔
1557
  const vector<ParentCell>& parent_cells() const { return parent_cells_; }
1558

1559
  // Data Members
1560
  vector<ParentCell> parent_cells_;
1561
  std::unordered_map<int32_t, std::unordered_set<ParentCell, ParentCellHash>>
1562
    visited_cells_;
1563
};
1564

1565
vector<ParentCell> Cell::find_parent_cells(
×
1566
  int32_t instance, const Position& r) const
1567
{
1568

1569
  // create a temporary particle
1570
  GeometryState dummy_particle {};
×
1571
  dummy_particle.r() = r;
×
1572
  dummy_particle.u() = {0., 0., 1.};
×
1573

1574
  return find_parent_cells(instance, dummy_particle);
×
1575
}
×
1576

1577
vector<ParentCell> Cell::find_parent_cells(
×
1578
  int32_t instance, GeometryState& p) const
1579
{
1580
  // look up the particle's location
1581
  exhaustive_find_cell(p);
×
1582
  const auto& coords = p.coord();
×
1583

1584
  // build a parent cell stack from the particle coordinates
1585
  ParentCellStack stack;
×
1586
  bool cell_found = false;
×
1587
  for (auto it = coords.begin(); it != coords.end(); it++) {
×
1588
    const auto& coord = *it;
×
1589
    const auto& cell = model::cells[coord.cell()];
×
1590
    // if the cell at this level matches the current cell, stop adding to the
1591
    // stack
1592
    if (coord.cell() == model::cell_map[this->id_]) {
×
1593
      cell_found = true;
1594
      break;
1595
    }
1596

1597
    // if filled with a lattice, get the lattice index from the next
1598
    // level in the coordinates to push to the stack
1599
    int lattice_idx = C_NONE;
×
1600
    if (cell->type_ == Fill::LATTICE) {
×
1601
      const auto& next_coord = *(it + 1);
×
1602
      lattice_idx = model::lattices[next_coord.lattice()]->get_flat_index(
×
1603
        next_coord.lattice_index());
1604
    }
1605
    stack.push(coord.universe(), {coord.cell(), lattice_idx});
×
1606
  }
1607

1608
  // if this loop finished because the cell was found and
1609
  // the instance matches the one requested in the call
1610
  // we have the correct path and can return the stack
1611
  if (cell_found &&
×
1612
      stack.compute_instance(this->distribcell_index_) == instance) {
×
1613
    return stack.parent_cells();
×
1614
  }
1615

1616
  // fall back on an exhaustive search for the cell's parents
1617
  return exhaustive_find_parent_cells(instance);
×
1618
}
×
1619

1620
vector<ParentCell> Cell::exhaustive_find_parent_cells(int32_t instance) const
136✔
1621
{
1622
  ParentCellStack stack;
136✔
1623
  // start with this cell's universe
1624
  int32_t prev_univ_idx;
136✔
1625
  int32_t univ_idx = this->universe_;
136✔
1626

1627
  while (true) {
211✔
1628
    const auto& univ = model::universes[univ_idx];
211✔
1629
    prev_univ_idx = univ_idx;
211✔
1630

1631
    // search for a cell that is filled w/ this universe
1632
    for (const auto& cell : model::cells) {
1,429✔
1633
      // if this is a material-filled cell, move on
1634
      if (cell->type_ == Fill::MATERIAL)
1,323✔
1635
        continue;
782✔
1636

1637
      if (cell->type_ == Fill::UNIVERSE) {
541✔
1638
        // if this is in the set of cells previously visited for this universe,
1639
        // move on
1640
        if (stack.visited(univ_idx, {model::cell_map[cell->id_], C_NONE}))
316!
1641
          continue;
×
1642

1643
        // if this cell contains the universe we're searching for, add it to the
1644
        // stack
1645
        if (cell->fill_ == univ_idx) {
316!
1646
          stack.push(univ_idx, {model::cell_map[cell->id_], C_NONE});
×
1647
          univ_idx = cell->universe_;
×
1648
        }
1649
      } else if (cell->type_ == Fill::LATTICE) {
225!
1650
        // retrieve the lattice and lattice universes
1651
        const auto& lattice = model::lattices[cell->fill_];
225✔
1652
        const auto& lattice_univs = lattice->universes_;
225✔
1653

1654
        // start search for universe
1655
        auto lat_it = lattice_univs.begin();
225✔
1656
        while (true) {
495✔
1657
          // find the next lattice cell with this universe
1658
          lat_it = std::find(lat_it, lattice_univs.end(), univ_idx);
360✔
1659
          if (lat_it == lattice_univs.end())
360✔
1660
            break;
1661

1662
          int lattice_idx = lat_it - lattice_univs.begin();
240✔
1663

1664
          // move iterator forward one to avoid finding the same entry
1665
          lat_it++;
240✔
1666
          if (stack.visited(
480✔
1667
                univ_idx, {model::cell_map[cell->id_], lattice_idx}))
240✔
1668
            continue;
135✔
1669

1670
          // add this cell and lattice index to the stack and exit loop
1671
          stack.push(univ_idx, {model::cell_map[cell->id_], lattice_idx});
105✔
1672
          univ_idx = cell->universe_;
105✔
1673
          break;
105✔
1674
        }
135✔
1675
      }
1676
      // if we've updated the universe, break
1677
      if (prev_univ_idx != univ_idx)
541✔
1678
        break;
1679
    } // end cell loop search for universe
1680

1681
    // if we're at the top of the geometry and the instance matches, we're done
1682
    if (univ_idx == model::root_universe &&
253!
1683
        stack.compute_instance(this->distribcell_index_) == instance)
211✔
1684
      break;
1685

1686
    // if there is no match on the original cell's universe, report an error
1687
    if (univ_idx == this->universe_) {
75!
1688
      fatal_error(
×
1689
        fmt::format("Could not find the parent cells for cell {}, instance {}.",
×
1690
          this->id_, instance));
×
1691
    }
1692

1693
    // if we don't find a suitable update, adjust the stack and continue
1694
    if (univ_idx == model::root_universe || univ_idx == prev_univ_idx) {
75!
1695
      stack.pop();
75✔
1696
      univ_idx = stack.empty() ? this->universe_ : stack.current_univ();
75!
1697
    }
1698

1699
  } // end while
1700

1701
  // reverse the stack so the highest cell comes first
1702
  std::reverse(stack.parent_cells().begin(), stack.parent_cells().end());
136✔
1703
  return stack.parent_cells();
272✔
1704
}
136✔
1705

1706
std::unordered_map<int32_t, vector<int32_t>> Cell::get_contained_cells(
181✔
1707
  int32_t instance, Position* hint) const
1708
{
1709
  std::unordered_map<int32_t, vector<int32_t>> contained_cells;
181✔
1710

1711
  // if this is a material-filled cell it has no contained cells
1712
  if (this->type_ == Fill::MATERIAL)
181✔
1713
    return contained_cells;
1714

1715
  // find the pathway through the geometry to this cell
1716
  vector<ParentCell> parent_cells;
136!
1717

1718
  // if a positional hint is provided, attempt to do a fast lookup
1719
  // of the parent cells
1720
  parent_cells = hint ? find_parent_cells(instance, *hint)
136!
1721
                      : exhaustive_find_parent_cells(instance);
136✔
1722

1723
  // if this cell is filled w/ a material, it contains no other cells
1724
  if (type_ != Fill::MATERIAL) {
136!
1725
    this->get_contained_cells_inner(contained_cells, parent_cells);
136✔
1726
  }
1727

1728
  return contained_cells;
136✔
1729
}
181✔
1730

1731
//! Get all cells within this cell
1732
void Cell::get_contained_cells_inner(
134,178✔
1733
  std::unordered_map<int32_t, vector<int32_t>>& contained_cells,
1734
  vector<ParentCell>& parent_cells) const
1735
{
1736

1737
  // filled by material, determine instance based on parent cells
1738
  if (type_ == Fill::MATERIAL) {
134,178✔
1739
    int instance = 0;
133,532✔
1740
    if (this->distribcell_index_ >= 0) {
133,532!
1741
      for (auto& parent_cell : parent_cells) {
400,594✔
1742
        auto& cell = model::cells[parent_cell.cell_index];
267,062✔
1743
        if (cell->type_ == Fill::UNIVERSE) {
267,062✔
1744
          instance += cell->offset_[distribcell_index_];
132,032✔
1745
        } else if (cell->type_ == Fill::LATTICE) {
135,030!
1746
          auto& lattice = model::lattices[cell->fill_];
135,030✔
1747
          instance += lattice->offset(
135,030✔
1748
            this->distribcell_index_, parent_cell.lattice_index);
135,030✔
1749
        }
1750
      }
1751
    }
1752
    // add entry to contained cells
1753
    contained_cells[model::cell_map[id_]].push_back(instance);
133,532✔
1754
    // filled with universe, add the containing cell to the parent cells
1755
    // and recurse
1756
  } else if (type_ == Fill::UNIVERSE) {
646✔
1757
    parent_cells.push_back({model::cell_map[id_], -1});
526✔
1758
    auto& univ = model::universes[fill_];
526✔
1759
    for (auto cell_index : univ->cells_) {
3,033✔
1760
      auto& cell = model::cells[cell_index];
2,507✔
1761
      cell->get_contained_cells_inner(contained_cells, parent_cells);
2,507✔
1762
    }
1763
    parent_cells.pop_back();
526✔
1764
    // filled with a lattice, visit each universe in the lattice
1765
    // with a recursive call to collect the cell instances
1766
  } else if (type_ == Fill::LATTICE) {
120!
1767
    auto& lattice = model::lattices[fill_];
120✔
1768
    for (auto i = lattice->begin(); i != lattice->end(); ++i) {
131,340✔
1769
      auto& univ = model::universes[*i];
131,220✔
1770
      parent_cells.push_back({model::cell_map[id_], i.indx_});
131,220✔
1771
      for (auto cell_index : univ->cells_) {
262,755✔
1772
        auto& cell = model::cells[cell_index];
131,535✔
1773
        cell->get_contained_cells_inner(contained_cells, parent_cells);
131,535✔
1774
      }
1775
      parent_cells.pop_back();
131,220✔
1776
    }
1777
  }
1778
}
134,178✔
1779

1780
//! Return the index in the cells array of a cell with a given ID
1781
extern "C" int openmc_get_cell_index(int32_t id, int32_t* index)
1,021✔
1782
{
1783
  auto it = model::cell_map.find(id);
1,021✔
1784
  if (it != model::cell_map.end()) {
1,021✔
1785
    *index = it->second;
1,010✔
1786
    return 0;
1,010✔
1787
  } else {
1788
    set_errmsg("No cell exists with ID=" + std::to_string(id) + ".");
11✔
1789
    return OPENMC_E_INVALID_ID;
11✔
1790
  }
1791
}
1792

1793
//! Return the ID of a cell
1794
extern "C" int openmc_cell_get_id(int32_t index, int32_t* id)
602,407✔
1795
{
1796
  if (index >= 0 && index < model::cells.size()) {
602,407!
1797
    *id = model::cells[index]->id_;
602,407✔
1798
    return 0;
602,407✔
1799
  } else {
1800
    set_errmsg("Index in cells array is out of bounds.");
×
1801
    return OPENMC_E_OUT_OF_BOUNDS;
×
1802
  }
1803
}
1804

1805
//! Set the ID of a cell
1806
extern "C" int openmc_cell_set_id(int32_t index, int32_t id)
22✔
1807
{
1808
  if (index >= 0 && index < model::cells.size()) {
22!
1809
    model::cells[index]->id_ = id;
22✔
1810
    model::cell_map[id] = index;
22✔
1811
    return 0;
22✔
1812
  } else {
1813
    set_errmsg("Index in cells array is out of bounds.");
×
1814
    return OPENMC_E_OUT_OF_BOUNDS;
×
1815
  }
1816
}
1817

1818
//! Return the translation vector of a cell
1819
extern "C" int openmc_cell_get_translation(int32_t index, double xyz[])
55✔
1820
{
1821
  if (index >= 0 && index < model::cells.size()) {
55!
1822
    auto& cell = model::cells[index];
55✔
1823
    xyz[0] = cell->translation_.x;
55✔
1824
    xyz[1] = cell->translation_.y;
55✔
1825
    xyz[2] = cell->translation_.z;
55✔
1826
    return 0;
55✔
1827
  } else {
1828
    set_errmsg("Index in cells array is out of bounds.");
×
1829
    return OPENMC_E_OUT_OF_BOUNDS;
×
1830
  }
1831
}
1832

1833
//! Set the translation vector of a cell
1834
extern "C" int openmc_cell_set_translation(int32_t index, const double xyz[])
55✔
1835
{
1836
  if (index >= 0 && index < model::cells.size()) {
55!
1837
    if (model::cells[index]->fill_ == C_NONE) {
55✔
1838
      set_errmsg(fmt::format("Cannot apply a translation to cell {}"
11✔
1839
                             " because it is not filled with another universe",
1840
        index));
1841
      return OPENMC_E_GEOMETRY;
11✔
1842
    }
1843
    model::cells[index]->translation_ = Position(xyz);
44✔
1844
    return 0;
44✔
1845
  } else {
1846
    set_errmsg("Index in cells array is out of bounds.");
×
1847
    return OPENMC_E_OUT_OF_BOUNDS;
×
1848
  }
1849
}
1850

1851
//! Return the rotation matrix of a cell
1852
extern "C" int openmc_cell_get_rotation(int32_t index, double rot[], size_t* n)
55✔
1853
{
1854
  if (index >= 0 && index < model::cells.size()) {
55!
1855
    auto& cell = model::cells[index];
55✔
1856
    *n = cell->rotation_.size();
55✔
1857
    std::memcpy(rot, cell->rotation_.data(), *n * sizeof(cell->rotation_[0]));
55✔
1858
    return 0;
55✔
1859
  } else {
1860
    set_errmsg("Index in cells array is out of bounds.");
×
1861
    return OPENMC_E_OUT_OF_BOUNDS;
×
1862
  }
1863
}
1864

1865
//! Set the flattened rotation matrix of a cell
1866
extern "C" int openmc_cell_set_rotation(
66✔
1867
  int32_t index, const double rot[], size_t rot_len)
1868
{
1869
  if (index >= 0 && index < model::cells.size()) {
66!
1870
    if (model::cells[index]->fill_ == C_NONE) {
66✔
1871
      set_errmsg(fmt::format("Cannot apply a rotation to cell {}"
11✔
1872
                             " because it is not filled with another universe",
1873
        index));
1874
      return OPENMC_E_GEOMETRY;
11✔
1875
    }
1876
    std::vector<double> vec_rot(rot, rot + rot_len);
55✔
1877
    model::cells[index]->set_rotation(vec_rot);
55✔
1878
    return 0;
55✔
1879
  } else {
66✔
1880
    set_errmsg("Index in cells array is out of bounds.");
×
1881
    return OPENMC_E_OUT_OF_BOUNDS;
×
1882
  }
1883
}
1884

1885
//! Get the number of instances of the requested cell
1886
extern "C" int openmc_cell_get_num_instances(
77✔
1887
  int32_t index, int32_t* num_instances)
1888
{
1889
  if (index < 0 || index >= model::cells.size()) {
77!
1890
    set_errmsg("Index in cells array is out of bounds.");
×
1891
    return OPENMC_E_OUT_OF_BOUNDS;
×
1892
  }
1893
  *num_instances = model::cells[index]->n_instances();
77✔
1894
  return 0;
77✔
1895
}
1896

1897
//! Extend the cells array by n elements
1898
extern "C" int openmc_extend_cells(
22✔
1899
  int32_t n, int32_t* index_start, int32_t* index_end)
1900
{
1901
  if (index_start)
22!
1902
    *index_start = model::cells.size();
22✔
1903
  if (index_end)
22!
1904
    *index_end = model::cells.size() + n - 1;
×
1905
  for (int32_t i = 0; i < n; i++) {
44✔
1906
    model::cells.push_back(make_unique<CSGCell>());
22✔
1907
  }
1908
  return 0;
22✔
1909
}
1910

1911
extern "C" int cells_size()
99✔
1912
{
1913
  return model::cells.size();
99✔
1914
}
1915

1916
} // 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