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

openmc-dev / openmc / 20470713200

23 Dec 2025 08:18PM UTC coverage: 81.94% (-0.2%) from 82.144%
20470713200

push

github

web-flow
Refactor get_energy_index to prevent repetition (#3686)

17063 of 23868 branches covered (71.49%)

Branch coverage included in aggregate %.

11 of 11 new or added lines in 4 files covered. (100.0%)

4 existing lines in 2 files now uncovered.

55178 of 64295 relevant lines covered (85.82%)

43466998.05 hits per line

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

65.95
/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
4,240✔
44
{
45
  return model::universes[universe_]->n_instances_;
4,240✔
46
}
47

48
void Cell::set_rotation(const vector<double>& rot)
394✔
49
{
50
  if (fill_ == C_NONE) {
394!
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) {
394!
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();
394✔
62
  rotation_.reserve(rot.size() == 9 ? 9 : 12);
394!
63
  if (rot.size() == 3) {
394!
64
    double phi = -rot[0] * PI / 180.0;
394✔
65
    double theta = -rot[1] * PI / 180.0;
394✔
66
    double psi = -rot[2] * PI / 180.0;
394✔
67
    rotation_.push_back(std::cos(theta) * std::cos(psi));
394!
68
    rotation_.push_back(-std::cos(phi) * std::sin(psi) +
×
69
                        std::sin(phi) * std::sin(theta) * std::cos(psi));
394!
70
    rotation_.push_back(std::sin(phi) * std::sin(psi) +
×
71
                        std::cos(phi) * std::sin(theta) * std::cos(psi));
394!
72
    rotation_.push_back(std::cos(theta) * std::sin(psi));
394!
73
    rotation_.push_back(std::cos(phi) * std::cos(psi) +
×
74
                        std::sin(phi) * std::sin(theta) * std::sin(psi));
394!
75
    rotation_.push_back(-std::sin(phi) * std::cos(psi) +
×
76
                        std::cos(phi) * std::sin(theta) * std::sin(psi));
394!
77
    rotation_.push_back(-std::sin(theta));
394!
78
    rotation_.push_back(std::sin(phi) * std::cos(theta));
394!
79
    rotation_.push_back(std::cos(phi) * std::cos(theta));
394!
80

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

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

96
  if (instance >= 0) {
80✔
97
    double sqrtkT = sqrtkT_.size() == 1 ? sqrtkT_.at(0) : sqrtkT_.at(instance);
10!
98
    return sqrtkT * sqrtkT / K_BOLTZMANN;
10✔
99
  } else {
100
    return sqrtkT_[0] * sqrtkT_[0] / K_BOLTZMANN;
70✔
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);
2,147,483,647✔
109
  } else {
110
    return density_mult_[0];
70✔
111
  }
112
}
113

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

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

123
void Cell::set_temperature(double T, int32_t instance, bool set_contained)
438✔
124
{
125
  if (settings::temperature_method == TemperatureMethod::INTERPOLATION) {
438!
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) {
438✔
140
    if (instance >= 0) {
410✔
141
      // If temperature vector is not big enough, resize it first
142
      if (sqrtkT_.size() != n_instances())
340✔
143
        sqrtkT_.resize(n_instances(), sqrtkT_[0]);
42✔
144

145
      // Set temperature for the corresponding instance
146
      sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T);
340✔
147
    } else {
148
      // Set temperature for all instances
149
      for (auto& T_ : sqrtkT_) {
140✔
150
        T_ = std::sqrt(K_BOLTZMANN * T);
70✔
151
      }
152
    }
153
  } else {
154
    if (!set_contained) {
28!
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);
28!
162
    for (const auto& entry : contained_cells) {
112✔
163
      auto& cell = model::cells[entry.first];
84✔
164
      assert(cell->type_ == Fill::MATERIAL);
66!
165
      auto& instances = entry.second;
84✔
166
      for (auto instance : instances) {
294✔
167
        cell->set_temperature(T, instance);
210!
168
      }
169
    }
170
  }
28✔
171
}
438✔
172

173
void Cell::set_density(double density, int32_t instance, bool set_contained)
256✔
174
{
175
  if (type_ != Fill::MATERIAL && !set_contained) {
256!
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) {
256✔
183
    const int32_t mat_index = material(instance);
242✔
184
    if (mat_index == MATERIAL_VOID)
242!
185
      return;
×
186

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

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

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

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

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

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

234
  close_group(cell_group);
140!
235
}
140✔
236

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

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

245
  // Ensure number of temperatures makes sense
246
  auto n_temps = temps.size();
160✔
247
  if (n_temps > 1 && n_temps != n_instances()) {
160!
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();
160✔
255
  sqrtkT_.resize(temps.size());
160!
256
  for (int64_t i = 0; i < temps.size(); ++i) {
280✔
257
    this->set_temperature(temps[i], i);
120!
258
  }
259

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

265
    // Ensure number of densities makes sense
266
    auto n_density = density.size();
120✔
267
    if (n_density > 1 && n_density != n_instances()) {
120!
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) {
240✔
275
      this->set_density(density[i], i);
120!
276
    }
277
  }
120✔
278

279
  close_group(cell_group);
160!
280
}
160✔
281

282
void Cell::to_hdf5(hid_t cell_group) const
22,695✔
283
{
284

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

288
  if (!name_.empty()) {
22,695✔
289
    write_string(group, "name", name_, false);
4,672✔
290
  }
291

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

294
  to_hdf5_inner(group);
22,695✔
295

296
  // Write fill information.
297
  if (type_ == Fill::MATERIAL) {
22,695✔
298
    write_dataset(group, "fill_type", "material");
18,155!
299
    std::vector<int32_t> mat_ids;
18,155✔
300
    for (auto i_mat : material_) {
37,454✔
301
      if (i_mat != MATERIAL_VOID) {
19,299✔
302
        mat_ids.push_back(model::materials[i_mat]->id_);
13,438!
303
      } else {
304
        mat_ids.push_back(MATERIAL_VOID);
5,861!
305
      }
306
    }
307
    if (mat_ids.size() == 1) {
18,155✔
308
      write_dataset(group, "material", mat_ids[0]);
17,987!
309
    } else {
310
      write_dataset(group, "material", mat_ids);
168!
311
    }
312

313
    std::vector<double> temps;
18,155✔
314
    for (auto sqrtkT_val : sqrtkT_)
37,574✔
315
      temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN);
19,419!
316
    write_dataset(group, "temperature", temps);
18,155!
317

318
    write_dataset(group, "density_mult", density_mult_);
18,155!
319

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

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

340
  close_group(group);
22,695✔
341
}
22,695✔
342

343
//==============================================================================
344
// CSGCell implementation
345
//==============================================================================
346

347
CSGCell::CSGCell(pugi::xml_node cell_node)
27,994✔
348
{
349
  if (check_for_node(cell_node, "id")) {
27,994!
350
    id_ = std::stoi(get_node_value(cell_node, "id"));
27,994!
351
  } else {
352
    fatal_error("Must specify id of cell in geometry XML file.");
×
353
  }
354

355
  if (check_for_node(cell_node, "name")) {
27,994!
356
    name_ = get_node_value(cell_node, "name");
6,287!
357
  }
358

359
  if (check_for_node(cell_node, "universe")) {
27,994!
360
    universe_ = std::stoi(get_node_value(cell_node, "universe"));
26,936!
361
  } else {
362
    universe_ = 0;
1,058✔
363
  }
364

365
  // Make sure that either material or fill was specified, but not both.
366
  bool fill_present = check_for_node(cell_node, "fill");
27,994!
367
  bool material_present = check_for_node(cell_node, "material");
27,994!
368
  if (!(fill_present || material_present)) {
27,994!
369
    fatal_error(
×
370
      fmt::format("Neither material nor fill was specified for cell {}", id_));
×
371
  }
372
  if (fill_present && material_present) {
27,994!
373
    fatal_error(fmt::format("Cell {} has both a material and a fill specified; "
×
374
                            "only one can be specified per cell",
375
      id_));
×
376
  }
377

378
  if (fill_present) {
27,994✔
379
    fill_ = std::stoi(get_node_value(cell_node, "fill"));
6,128!
380
    if (fill_ == universe_) {
6,128!
381
      fatal_error(fmt::format("Cell {} is filled with the same universe that "
×
382
                              "it is contained in.",
383
        id_));
×
384
    }
385
  } else {
386
    fill_ = C_NONE;
21,866✔
387
  }
388

389
  // Read the material element.  There can be zero materials (filled with a
390
  // universe), more than one material (distribmats), and some materials may
391
  // be "void".
392
  if (material_present) {
27,994✔
393
    vector<std::string> mats {
394
      get_node_array<std::string>(cell_node, "material", true)};
21,866!
395
    if (mats.size() > 0) {
21,866!
396
      material_.reserve(mats.size());
21,866!
397
      for (std::string mat : mats) {
44,912!
398
        if (mat.compare("void") == 0) {
23,046✔
399
          material_.push_back(MATERIAL_VOID);
6,233!
400
        } else {
401
          material_.push_back(std::stoi(mat));
16,813!
402
        }
403
      }
23,046✔
404
    } else {
405
      fatal_error(fmt::format(
×
406
        "An empty material element was specified for cell {}", id_));
×
407
    }
408
  }
21,866✔
409

410
  // Read the temperature element which may be distributed like materials.
411
  if (check_for_node(cell_node, "temperature")) {
27,994!
412
    sqrtkT_ = get_node_array<double>(cell_node, "temperature");
276!
413
    sqrtkT_.shrink_to_fit();
276!
414

415
    // Make sure this is a material-filled cell.
416
    if (material_.size() == 0) {
276!
417
      fatal_error(fmt::format(
×
418
        "Cell {} was specified with a temperature but no material. Temperature"
419
        "specification is only valid for cells filled with a material.",
420
        id_));
×
421
    }
422

423
    // Make sure all temperatures are non-negative.
424
    for (auto T : sqrtkT_) {
594✔
425
      if (T < 0) {
318!
426
        fatal_error(fmt::format(
×
427
          "Cell {} was specified with a negative temperature", id_));
×
428
      }
429
    }
430

431
    // Convert to sqrt(k*T).
432
    for (auto& T : sqrtkT_) {
594✔
433
      T = std::sqrt(K_BOLTZMANN * T);
318✔
434
    }
435
  }
436

437
  // Read the density element which can be distributed similar to temperature.
438
  // These get assigned to the density multiplier, requiring a division by
439
  // the material density.
440
  // Note: calculating the actual density multiplier is deferred until materials
441
  // are finalized. density_mult_ contains the true density in the meantime.
442
  if (check_for_node(cell_node, "density")) {
27,994!
443
    density_mult_ = get_node_array<double>(cell_node, "density");
14!
444
    density_mult_.shrink_to_fit();
14!
445

446
    // Make sure this is a material-filled cell.
447
    if (material_.size() == 0) {
14!
448
      fatal_error(fmt::format(
×
449
        "Cell {} was specified with a density but no material. Density"
450
        "specification is only valid for cells filled with a material.",
451
        id_));
×
452
    }
453

454
    // Make sure this is a non-void material.
455
    for (auto mat_id : material_) {
28✔
456
      if (mat_id == MATERIAL_VOID) {
14!
457
        fatal_error(fmt::format(
×
458
          "Cell {} was specified with a density, but contains a void "
459
          "material. Density specification is only valid for cells "
460
          "filled with a non-void material.",
461
          id_));
×
462
      }
463
    }
464

465
    // Make sure all densities are non-negative and greater than zero.
466
    for (auto rho : density_mult_) {
70✔
467
      if (rho <= 0) {
56!
468
        fatal_error(fmt::format(
×
469
          "Cell {} was specified with a density less than or equal to zero",
470
          id_));
×
471
      }
472
    }
473
  }
474

475
  // Read the region specification.
476
  std::string region_spec;
27,994✔
477
  if (check_for_node(cell_node, "region")) {
27,994!
478
    region_spec = get_node_value(cell_node, "region");
21,456!
479
  }
480

481
  // Get a tokenized representation of the region specification and apply De
482
  // Morgans law
483
  Region region(region_spec, id_);
27,994!
484
  region_ = region;
27,994!
485

486
  // Read the translation vector.
487
  if (check_for_node(cell_node, "translation")) {
27,994!
488
    if (fill_ == C_NONE) {
2,386!
489
      fatal_error(fmt::format("Cannot apply a translation to cell {}"
×
490
                              " because it is not filled with another universe",
491
        id_));
×
492
    }
493

494
    auto xyz {get_node_array<double>(cell_node, "translation")};
2,386!
495
    if (xyz.size() != 3) {
2,386!
496
      fatal_error(
×
497
        fmt::format("Non-3D translation vector applied to cell {}", id_));
×
498
    }
499
    translation_ = xyz;
2,386✔
500
  }
2,386✔
501

502
  // Read the rotation transform.
503
  if (check_for_node(cell_node, "rotation")) {
27,994!
504
    auto rot {get_node_array<double>(cell_node, "rotation")};
364!
505
    set_rotation(rot);
364!
506
  }
364✔
507
}
27,994✔
508

509
//==============================================================================
510

511
void CSGCell::to_hdf5_inner(hid_t group_id) const
22,695✔
512
{
513
  write_string(group_id, "geom_type", "csg", false);
22,695!
514
  write_string(group_id, "region", region_.str(), false);
22,695!
515
}
22,695✔
516

517
//==============================================================================
518

519
vector<int32_t>::iterator CSGCell::find_left_parenthesis(
×
520
  vector<int32_t>::iterator start, const vector<int32_t>& infix)
521
{
522
  // start search at zero
523
  int parenthesis_level = 0;
×
524
  auto it = start;
×
525
  while (it != infix.begin()) {
×
526
    // look at two tokens at a time
527
    int32_t one = *it;
×
528
    int32_t two = *(it - 1);
×
529

530
    // decrement parenthesis level if there are two adjacent surfaces
531
    if (one < OP_UNION && two < OP_UNION) {
×
532
      parenthesis_level--;
×
533
      // increment if there are two adjacent operators
534
    } else if (one >= OP_UNION && two >= OP_UNION) {
×
535
      parenthesis_level++;
×
536
    }
537

538
    // if the level gets to zero, return the position
539
    if (parenthesis_level == 0) {
×
540
      // move the iterator back one before leaving the loop
541
      // so that all tokens in the parenthesis block are included
542
      it--;
×
543
      break;
×
544
    }
545

546
    // continue loop, one token at a time
547
    it--;
×
548
  }
549
  return it;
×
550
}
551

552
//==============================================================================
553
// Region implementation
554
//==============================================================================
555

556
Region::Region(std::string region_spec, int32_t cell_id)
27,994✔
557
{
558
  // Check if region_spec is not empty.
559
  if (!region_spec.empty()) {
27,994✔
560
    // Parse all halfspaces and operators except for intersection (whitespace).
561
    for (int i = 0; i < region_spec.size();) {
124,746✔
562
      if (region_spec[i] == '(') {
103,290!
563
        expression_.push_back(OP_LEFT_PAREN);
1,027!
564
        i++;
1,027✔
565

566
      } else if (region_spec[i] == ')') {
102,263!
567
        expression_.push_back(OP_RIGHT_PAREN);
1,027!
568
        i++;
1,027✔
569

570
      } else if (region_spec[i] == '|') {
101,236!
571
        expression_.push_back(OP_UNION);
3,190!
572
        i++;
3,190✔
573

574
      } else if (region_spec[i] == '~') {
98,046!
575
        expression_.push_back(OP_COMPLEMENT);
28!
576
        i++;
28✔
577

578
      } else if (region_spec[i] == '-' || region_spec[i] == '+' ||
165,274!
579
                 std::isdigit(region_spec[i])) {
67,256!
580
        // This is the start of a halfspace specification.  Iterate j until we
581
        // find the end, then push-back everything between i and j.
582
        int j = i + 1;
57,446✔
583
        while (j < region_spec.size() && std::isdigit(region_spec[j])) {
117,260!
584
          j++;
59,814✔
585
        }
586
        expression_.push_back(std::stoi(region_spec.substr(i, j - i)));
57,446!
587
        i = j;
57,446✔
588

589
      } else if (std::isspace(region_spec[i])) {
40,572!
590
        i++;
40,572✔
591

592
      } else {
593
        auto err_msg =
594
          fmt::format("Region specification contains invalid character, \"{}\"",
595
            region_spec[i]);
×
596
        fatal_error(err_msg);
×
597
      }
×
598
    }
599

600
    // Add in intersection operators where a missing operator is needed.
601
    int i = 0;
21,456✔
602
    while (i < expression_.size() - 1) {
95,518✔
603
      bool left_compat {
604
        (expression_[i] < OP_UNION) || (expression_[i] == OP_RIGHT_PAREN)};
74,062✔
605
      bool right_compat {(expression_[i + 1] < OP_UNION) ||
74,062✔
606
                         (expression_[i + 1] == OP_LEFT_PAREN) ||
78,335✔
607
                         (expression_[i + 1] == OP_COMPLEMENT)};
4,273✔
608
      if (left_compat && right_compat) {
74,062✔
609
        expression_.insert(expression_.begin() + i + 1, OP_INTERSECTION);
32,800!
610
      }
611
      i++;
74,062✔
612
    }
613

614
    // Remove complement operators using DeMorgan's laws
615
    auto it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT);
21,456!
616
    while (it != expression_.end()) {
21,484✔
617
      // Erase complement
618
      expression_.erase(it);
28!
619

620
      // Define stop given left parenthesis or not
621
      auto stop = it;
28✔
622
      if (*it == OP_LEFT_PAREN) {
28!
623
        int depth = 1;
28✔
624
        do {
625
          stop++;
224✔
626
          if (*stop > OP_COMPLEMENT) {
224✔
627
            if (*stop == OP_RIGHT_PAREN) {
28!
628
              depth--;
28✔
629
            } else {
630
              depth++;
×
631
            }
632
          }
633
        } while (depth > 0);
224✔
634
        it++;
28✔
635
      }
636

637
      // apply DeMorgan's law to any surfaces/operators between these
638
      // positions in the RPN
639
      apply_demorgan(it, stop);
28!
640
      // update iterator position
641
      it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT);
28!
642
    }
643

644
    // Convert user IDs to surface indices.
645
    for (auto& r : expression_) {
116,946✔
646
      if (r < OP_UNION) {
95,490✔
647
        const auto& it {model::surface_map.find(abs(r))};
57,446!
648
        if (it == model::surface_map.end()) {
57,446!
649
          throw std::runtime_error {
×
650
            "Invalid surface ID " + std::to_string(abs(r)) +
×
651
            " specified in region for cell " + std::to_string(cell_id) + "."};
×
652
        }
653
        r = (r > 0) ? it->second + 1 : -(it->second + 1);
57,446✔
654
      }
655
    }
656

657
    // Check if this is a simple cell.
658
    simple_ = true;
21,456✔
659
    for (int32_t token : expression_) {
106,312✔
660
      if (token == OP_UNION) {
85,818✔
661
        simple_ = false;
962✔
662
        // Ensure intersections have precedence over unions
663
        add_precedence();
962!
664
        break;
962✔
665
      }
666
    }
667

668
    // If this cell is simple, remove all the superfluous operator tokens.
669
    if (simple_) {
21,456✔
670
      for (auto it = expression_.begin(); it != expression_.end(); it++) {
99,936✔
671
        if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) {
79,442!
672
          expression_.erase(it);
29,474!
673
          it--;
29,474✔
674
        }
675
      }
676
    }
677
    expression_.shrink_to_fit();
21,456!
678

679
  } else {
680
    simple_ = true;
6,538✔
681
  }
682
}
27,994✔
683

684
//==============================================================================
685

686
void Region::apply_demorgan(
196✔
687
  vector<int32_t>::iterator start, vector<int32_t>::iterator stop)
688
{
689
  do {
690
    if (*start < OP_UNION) {
196✔
691
      *start *= -1;
112✔
692
    } else if (*start == OP_UNION) {
84!
693
      *start = OP_INTERSECTION;
×
694
    } else if (*start == OP_INTERSECTION) {
84!
695
      *start = OP_UNION;
84✔
696
    }
697
    start++;
196✔
698
  } while (start < stop);
196✔
699
}
28✔
700

701
//==============================================================================
702
//! Add precedence for infix regions so intersections have higher
703
//! precedence than unions using parentheses.
704
//==============================================================================
705

706
int64_t Region::add_parentheses(int64_t start)
28✔
707
{
708
  int32_t start_token = expression_[start];
28✔
709
  // Add left parenthesis and set new position to be after parenthesis
710
  if (start_token == OP_UNION) {
28!
711
    start += 2;
×
712
  }
713
  expression_.insert(expression_.begin() + start - 1, OP_LEFT_PAREN);
28!
714

715
  // Keep track of return iterator distance. If we don't encounter a left
716
  // parenthesis, we return an iterator corresponding to wherever the right
717
  // parenthesis is inserted. If a left parenthesis is encountered, an iterator
718
  // corresponding to the left parenthesis is returned. Also note that we keep
719
  // track of a *distance* instead of an iterator because the underlying memory
720
  // allocation may change.
721
  std::size_t return_it_dist = 0;
28✔
722

723
  // Add right parenthesis
724
  // While the start iterator is within the bounds of infix
725
  while (start + 1 < expression_.size()) {
196!
726
    start++;
196✔
727

728
    // If the current token is an operator and is different than the start token
729
    if (expression_[start] >= OP_UNION && expression_[start] != start_token) {
196✔
730
      // Skip wrapped regions but save iterator position to check precedence and
731
      // add right parenthesis, right parenthesis position depends on the
732
      // operator, when the operator is a union then do not include the operator
733
      // in the region, when the operator is an intersection then include the
734
      // operator and next surface
735
      if (expression_[start] == OP_LEFT_PAREN) {
28!
736
        return_it_dist = start;
×
737
        int depth = 1;
×
738
        do {
739
          start++;
×
740
          if (expression_[start] > OP_COMPLEMENT) {
×
741
            if (expression_[start] == OP_RIGHT_PAREN) {
×
742
              depth--;
×
743
            } else {
744
              depth++;
×
745
            }
746
          }
747
        } while (depth > 0);
×
748
      } else {
749
        if (start_token == OP_UNION) {
28!
750
          --start;
×
751
        }
752
        expression_.insert(expression_.begin() + start, OP_RIGHT_PAREN);
28!
753
        if (return_it_dist > 0) {
28!
754
          return return_it_dist;
×
755
        } else {
756
          return start - 1;
28✔
757
        }
758
      }
759
    }
760
  }
761
  // If we get here a right parenthesis hasn't been placed,
762
  // return iterator
763
  expression_.push_back(OP_RIGHT_PAREN);
×
764
  if (return_it_dist > 0) {
×
765
    return return_it_dist;
×
766
  } else {
767
    return start - 1;
×
768
  }
769
}
770

771
//==============================================================================
772

773
void Region::add_precedence()
962✔
774
{
775
  int32_t current_op = 0;
962✔
776
  std::size_t current_dist = 0;
962✔
777

778
  for (int64_t i = 0; i < expression_.size(); i++) {
16,982✔
779
    int32_t token = expression_[i];
16,020✔
780

781
    if (token == OP_UNION || token == OP_INTERSECTION) {
16,020✔
782
      if (current_op == 0) {
6,502✔
783
        // Set the current operator if is hasn't been set
784
        current_op = token;
2,073✔
785
        current_dist = i;
2,073✔
786
      } else if (token != current_op) {
4,429✔
787
        // If the current operator doesn't match the token, add parenthesis to
788
        // assert precedence
789
        if (current_op == OP_INTERSECTION) {
28✔
790
          i = add_parentheses(current_dist);
14✔
791
        } else {
792
          i = add_parentheses(i);
14✔
793
        }
794
        current_op = 0;
28✔
795
        current_dist = 0;
28✔
796
      }
797
    } else if (token > OP_COMPLEMENT) {
9,518✔
798
      // If the token is a parenthesis reset the current operator
799
      current_op = 0;
2,082✔
800
      current_dist = 0;
2,082✔
801
    }
802
  }
803
}
962✔
804

805
//==============================================================================
806
//! Convert infix region specification to Reverse Polish Notation (RPN)
807
//!
808
//! This function uses the shunting-yard algorithm.
809
//==============================================================================
810

811
vector<int32_t> Region::generate_postfix(int32_t cell_id) const
40✔
812
{
813
  vector<int32_t> rpn;
40✔
814
  vector<int32_t> stack;
40✔
815

816
  for (int32_t token : expression_) {
900✔
817
    if (token < OP_UNION) {
860✔
818
      // If token is not an operator, add it to output
819
      rpn.push_back(token);
360!
820
    } else if (token < OP_RIGHT_PAREN) {
500✔
821
      // Regular operators union, intersection, complement
822
      while (stack.size() > 0) {
510✔
823
        int32_t op = stack.back();
420✔
824

825
        if (op < OP_RIGHT_PAREN && ((token == OP_COMPLEMENT && token < op) ||
420!
826
                                     (token != OP_COMPLEMENT && token <= op))) {
190!
827
          // While there is an operator, op, on top of the stack, if the token
828
          // is left-associative and its precedence is less than or equal to
829
          // that of op or if the token is right-associative and its precedence
830
          // is less than that of op, move op to the output queue and push the
831
          // token on to the stack. Note that only complement is
832
          // right-associative.
833
          rpn.push_back(op);
190!
834
          stack.pop_back();
190✔
835
        } else {
836
          break;
837
        }
838
      }
839

840
      stack.push_back(token);
320!
841

842
    } else if (token == OP_LEFT_PAREN) {
180✔
843
      // If the token is a left parenthesis, push it onto the stack
844
      stack.push_back(token);
90!
845

846
    } else {
847
      // If the token is a right parenthesis, move operators from the stack to
848
      // the output queue until reaching the left parenthesis.
849
      for (auto it = stack.rbegin(); *it != OP_LEFT_PAREN; it++) {
180!
850
        // If we run out of operators without finding a left parenthesis, it
851
        // means there are mismatched parentheses.
852
        if (it == stack.rend()) {
90!
853
          fatal_error(fmt::format(
×
854
            "Mismatched parentheses in region specification for cell {}",
855
            cell_id));
856
        }
857
        rpn.push_back(stack.back());
90!
858
        stack.pop_back();
90✔
859
      }
860

861
      // Pop the left parenthesis.
862
      stack.pop_back();
90✔
863
    }
864
  }
865

866
  while (stack.size() > 0) {
80✔
867
    int32_t op = stack.back();
40✔
868

869
    // If the operator is a parenthesis it is mismatched.
870
    if (op >= OP_RIGHT_PAREN) {
40!
871
      fatal_error(fmt::format(
×
872
        "Mismatched parentheses in region specification for cell {}", cell_id));
873
    }
874

875
    rpn.push_back(stack.back());
40!
876
    stack.pop_back();
40✔
877
  }
878

879
  return rpn;
80✔
880
}
40✔
881

882
//==============================================================================
883

884
std::string Region::str() const
22,695✔
885
{
886
  std::stringstream region_spec {};
22,695!
887
  if (!expression_.empty()) {
22,695✔
888
    for (int32_t token : expression_) {
70,103✔
889
      if (token == OP_LEFT_PAREN) {
53,304✔
890
        region_spec << " (";
927!
891
      } else if (token == OP_RIGHT_PAREN) {
52,377✔
892
        region_spec << " )";
927!
893
      } else if (token == OP_COMPLEMENT) {
51,450!
894
        region_spec << " ~";
×
895
      } else if (token == OP_INTERSECTION) {
51,450✔
896
      } else if (token == OP_UNION) {
48,652✔
897
        region_spec << " |";
2,880!
898
      } else {
899
        // Note the off-by-one indexing
900
        auto surf_id = model::surfaces[abs(token) - 1]->id_;
45,772✔
901
        region_spec << " " << ((token > 0) ? surf_id : -surf_id);
45,772!
902
      }
903
    }
904
  }
905
  return region_spec.str();
45,390!
906
}
22,695✔
907

908
//==============================================================================
909

910
std::pair<double, int32_t> Region::distance(
2,147,483,647✔
911
  Position r, Direction u, int32_t on_surface) const
912
{
913
  double min_dist {INFTY};
2,147,483,647✔
914
  int32_t i_surf {std::numeric_limits<int32_t>::max()};
2,147,483,647✔
915

916
  for (int32_t token : expression_) {
2,147,483,647✔
917
    // Ignore this token if it corresponds to an operator rather than a region.
918
    if (token >= OP_UNION)
2,147,483,647✔
919
      continue;
124,365,151✔
920

921
    // Calculate the distance to this surface.
922
    // Note the off-by-one indexing
923
    bool coincident {std::abs(token) == std::abs(on_surface)};
2,147,483,647✔
924
    double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
2,147,483,647!
925

926
    // Check if this distance is the new minimum.
927
    if (d < min_dist) {
2,147,483,647✔
928
      if (min_dist - d >= FP_PRECISION * min_dist) {
2,147,483,647!
929
        min_dist = d;
2,147,483,647✔
930
        i_surf = -token;
2,147,483,647✔
931
      }
932
    }
933
  }
934

935
  return {min_dist, i_surf};
2,147,483,647✔
936
}
937

938
//==============================================================================
939

940
bool Region::contains(Position r, Direction u, int32_t on_surface) const
2,147,483,647✔
941
{
942
  if (simple_) {
2,147,483,647✔
943
    return contains_simple(r, u, on_surface);
2,147,483,647✔
944
  } else {
945
    return contains_complex(r, u, on_surface);
5,716,419✔
946
  }
947
}
948

949
//==============================================================================
950

951
bool Region::contains_simple(Position r, Direction u, int32_t on_surface) const
2,147,483,647✔
952
{
953
  for (int32_t token : expression_) {
2,147,483,647✔
954
    // Assume that no tokens are operators. Evaluate the sense of particle with
955
    // respect to the surface and see if the token matches the sense. If the
956
    // particle's surface attribute is set and matches the token, that
957
    // overrides the determination based on sense().
958
    if (token == on_surface) {
2,147,483,647✔
959
    } else if (-token == on_surface) {
2,147,483,647✔
960
      return false;
893,238,257✔
961
    } else {
962
      // Note the off-by-one indexing
963
      bool sense = model::surfaces[abs(token) - 1]->sense(r, u);
2,147,483,647!
964
      if (sense != (token > 0)) {
2,147,483,647✔
965
        return false;
759,718,930✔
966
      }
967
    }
968
  }
969
  return true;
2,147,483,647✔
970
}
971

972
//==============================================================================
973

974
bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const
5,716,419✔
975
{
976
  bool in_cell = true;
5,716,419✔
977
  int total_depth = 0;
5,716,419✔
978

979
  // For each token
980
  for (auto it = expression_.begin(); it != expression_.end(); it++) {
87,809,084✔
981
    int32_t token = *it;
83,740,288✔
982

983
    // If the token is a surface evaluate the sense
984
    // If the token is a union or intersection check to
985
    // short circuit
986
    if (token < OP_UNION) {
83,740,288✔
987
      if (token == on_surface) {
36,702,673✔
988
        in_cell = true;
2,822,830✔
989
      } else if (-token == on_surface) {
33,879,843✔
990
        in_cell = false;
572,093✔
991
      } else {
992
        // Note the off-by-one indexing
993
        bool sense = model::surfaces[abs(token) - 1]->sense(r, u);
33,307,750!
994
        in_cell = (sense == (token > 0));
33,307,750✔
995
      }
996
    } else if ((token == OP_UNION && in_cell == true) ||
47,037,615✔
997
               (token == OP_INTERSECTION && in_cell == false)) {
23,301,475✔
998
      // If the total depth is zero return
999
      if (total_depth == 0) {
5,303,925✔
1000
        return in_cell;
1,647,623✔
1001
      }
1002

1003
      total_depth--;
3,656,302✔
1004

1005
      // While the iterator is within the bounds of the vector
1006
      int depth = 1;
3,656,302✔
1007
      do {
1008
        // Get next token
1009
        it++;
22,694,402✔
1010
        int32_t next_token = *it;
22,694,402✔
1011

1012
        // If the token is an a parenthesis
1013
        if (next_token > OP_COMPLEMENT) {
22,694,402✔
1014
          // Adjust depth accordingly
1015
          if (next_token == OP_RIGHT_PAREN) {
3,829,452✔
1016
            depth--;
3,742,877✔
1017
          } else {
1018
            depth++;
86,575✔
1019
          }
1020
        }
1021
      } while (depth > 0);
22,694,402✔
1022
    } else if (token == OP_LEFT_PAREN) {
45,389,992✔
1023
      total_depth++;
7,201,869✔
1024
    } else if (token == OP_RIGHT_PAREN) {
34,531,821✔
1025
      total_depth--;
3,545,567✔
1026
    }
1027
  }
1028
  return in_cell;
4,068,796✔
1029
}
1030

1031
//==============================================================================
1032

1033
BoundingBox Region::bounding_box(int32_t cell_id) const
80✔
1034
{
1035
  if (simple_) {
80✔
1036
    return bounding_box_simple();
40✔
1037
  } else {
1038
    auto postfix = generate_postfix(cell_id);
40!
1039
    return bounding_box_complex(postfix);
40!
1040
  }
40✔
1041
}
1042

1043
//==============================================================================
1044

1045
BoundingBox Region::bounding_box_simple() const
40✔
1046
{
1047
  BoundingBox bbox;
40✔
1048
  for (int32_t token : expression_) {
160✔
1049
    bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0);
120!
1050
  }
1051
  return bbox;
40✔
1052
}
1053

1054
//==============================================================================
1055

1056
BoundingBox Region::bounding_box_complex(vector<int32_t> postfix) const
40✔
1057
{
1058
  vector<BoundingBox> stack(postfix.size());
40!
1059
  int i_stack = -1;
40✔
1060

1061
  for (auto& token : postfix) {
720✔
1062
    if (token == OP_UNION) {
680✔
1063
      stack[i_stack - 1] = stack[i_stack - 1] | stack[i_stack];
140!
1064
      i_stack--;
140✔
1065
    } else if (token == OP_INTERSECTION) {
540✔
1066
      stack[i_stack - 1] = stack[i_stack - 1] & stack[i_stack];
180!
1067
      i_stack--;
180✔
1068
    } else {
1069
      i_stack++;
360✔
1070
      stack[i_stack] = model::surfaces[abs(token) - 1]->bounding_box(token > 0);
360!
1071
    }
1072
  }
1073

1074
  assert(i_stack == 0);
32!
1075
  return stack.front();
80✔
1076
}
40✔
1077

1078
//==============================================================================
1079

1080
vector<int32_t> Region::surfaces() const
4,670✔
1081
{
1082
  if (simple_) {
4,670!
1083
    return expression_;
4,670!
1084
  }
1085

1086
  vector<int32_t> surfaces = expression_;
×
1087

1088
  auto it = std::find_if(surfaces.begin(), surfaces.end(),
×
1089
    [&](const auto& value) { return value >= OP_UNION; });
×
1090

1091
  while (it != surfaces.end()) {
×
1092
    surfaces.erase(it);
×
1093

1094
    it = std::find_if(surfaces.begin(), surfaces.end(),
×
1095
      [&](const auto& value) { return value >= OP_UNION; });
×
1096
  }
1097

1098
  return surfaces;
×
1099
}
×
1100

1101
//==============================================================================
1102
// Non-method functions
1103
//==============================================================================
1104

1105
void read_cells(pugi::xml_node node)
6,736✔
1106
{
1107
  // Count the number of cells.
1108
  int n_cells = 0;
6,736✔
1109
  for (pugi::xml_node cell_node : node.children("cell")) {
34,730!
1110
    n_cells++;
27,994!
1111
  }
1112

1113
  // Loop over XML cell elements and populate the array.
1114
  model::cells.reserve(n_cells);
6,736✔
1115
  for (pugi::xml_node cell_node : node.children("cell")) {
34,730!
1116
    model::cells.push_back(make_unique<CSGCell>(cell_node));
27,994!
1117
  }
1118

1119
  // Fill the cell map.
1120
  for (int i = 0; i < model::cells.size(); i++) {
34,730✔
1121
    int32_t id = model::cells[i]->id_;
27,994✔
1122
    auto search = model::cell_map.find(id);
27,994!
1123
    if (search == model::cell_map.end()) {
27,994!
1124
      model::cell_map[id] = i;
27,994!
1125
    } else {
1126
      fatal_error(
×
1127
        fmt::format("Two or more cells use the same unique ID: {}", id));
×
1128
    }
1129
  }
1130

1131
  read_dagmc_universes(node);
6,736✔
1132

1133
  populate_universes();
6,736✔
1134

1135
  // Allocate the cell overlap count if necessary.
1136
  if (settings::check_overlaps) {
6,736✔
1137
    model::overlap_check_count.resize(model::cells.size(), 0);
108!
1138
  }
1139

1140
  if (model::cells.size() == 0) {
6,736!
1141
    fatal_error("No cells were found in the geometry.xml file");
×
1142
  }
1143
}
6,736✔
1144

1145
void populate_universes()
6,736✔
1146
{
1147
  // Used to map universe index to the index of an implicit complement cell for
1148
  // DAGMC universes
1149
  std::unordered_map<int, int> implicit_comp_cells;
6,736✔
1150

1151
  // Populate the Universe vector and map.
1152
  for (int index_cell = 0; index_cell < model::cells.size(); index_cell++) {
34,730✔
1153
    int32_t uid = model::cells[index_cell]->universe_;
27,994✔
1154
    auto it = model::universe_map.find(uid);
27,994!
1155
    if (it == model::universe_map.end()) {
27,994✔
1156
      model::universes.push_back(make_unique<Universe>());
15,594!
1157
      model::universes.back()->id_ = uid;
15,594✔
1158
      model::universes.back()->cells_.push_back(index_cell);
15,594!
1159
      model::universe_map[uid] = model::universes.size() - 1;
15,594!
1160
    } else {
1161
#ifdef OPENMC_DAGMC_ENABLED
1162
      // Skip implicit complement cells for now
1163
      Universe* univ = model::universes[it->second].get();
1164
      DAGUniverse* dag_univ = dynamic_cast<DAGUniverse*>(univ);
×
1165
      if (dag_univ && (dag_univ->implicit_complement_idx() == index_cell)) {
×
1166
        implicit_comp_cells[it->second] = index_cell;
×
1167
        continue;
1168
      }
1169
#endif
1170

1171
      model::universes[it->second]->cells_.push_back(index_cell);
12,400!
1172
    }
1173
  }
1174

1175
  // Add DAGUniverse implicit complement cells last
1176
  for (const auto& it : implicit_comp_cells) {
6,736!
UNCOV
1177
    int index_univ = it.first;
×
UNCOV
1178
    int index_cell = it.second;
×
UNCOV
1179
    model::universes[index_univ]->cells_.push_back(index_cell);
×
1180
  }
1181

1182
  model::universes.shrink_to_fit();
6,736!
1183
}
6,736✔
1184

1185
//==============================================================================
1186
// C-API functions
1187
//==============================================================================
1188

1189
extern "C" int openmc_cell_get_fill(
50✔
1190
  int32_t index, int* type, int32_t** indices, int32_t* n)
1191
{
1192
  if (index >= 0 && index < model::cells.size()) {
50!
1193
    Cell& c {*model::cells[index]};
50✔
1194
    *type = static_cast<int>(c.type_);
50✔
1195
    if (c.type_ == Fill::MATERIAL) {
50!
1196
      *indices = c.material_.data();
50✔
1197
      *n = c.material_.size();
50✔
1198
    } else {
1199
      *indices = &c.fill_;
×
1200
      *n = 1;
×
1201
    }
1202
  } else {
1203
    set_errmsg("Index in cells array is out of bounds.");
×
1204
    return OPENMC_E_OUT_OF_BOUNDS;
×
1205
  }
1206
  return 0;
50✔
1207
}
1208

1209
extern "C" int openmc_cell_set_fill(
10✔
1210
  int32_t index, int type, int32_t n, const int32_t* indices)
1211
{
1212
  Fill filltype = static_cast<Fill>(type);
10✔
1213
  if (index >= 0 && index < model::cells.size()) {
10!
1214
    Cell& c {*model::cells[index]};
10✔
1215
    if (filltype == Fill::MATERIAL) {
10!
1216
      c.type_ = Fill::MATERIAL;
10✔
1217
      c.material_.clear();
10✔
1218
      for (int i = 0; i < n; i++) {
20✔
1219
        int i_mat = indices[i];
10✔
1220
        if (i_mat == MATERIAL_VOID) {
10!
1221
          c.material_.push_back(MATERIAL_VOID);
×
1222
        } else if (i_mat >= 0 && i_mat < model::materials.size()) {
10!
1223
          c.material_.push_back(i_mat);
10!
1224
        } else {
1225
          set_errmsg("Index in materials array is out of bounds.");
×
1226
          return OPENMC_E_OUT_OF_BOUNDS;
×
1227
        }
1228
      }
1229
      c.material_.shrink_to_fit();
10✔
1230
    } else if (filltype == Fill::UNIVERSE) {
×
1231
      c.type_ = Fill::UNIVERSE;
×
1232
    } else {
1233
      c.type_ = Fill::LATTICE;
×
1234
    }
1235
  } else {
1236
    set_errmsg("Index in cells array is out of bounds.");
×
1237
    return OPENMC_E_OUT_OF_BOUNDS;
×
1238
  }
1239
  return 0;
10✔
1240
}
1241

1242
extern "C" int openmc_cell_set_temperature(
80✔
1243
  int32_t index, double T, const int32_t* instance, bool set_contained)
1244
{
1245
  if (index < 0 || index >= model::cells.size()) {
80!
1246
    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
×
1247
    return OPENMC_E_OUT_OF_BOUNDS;
×
1248
  }
1249

1250
  int32_t instance_index = instance ? *instance : -1;
80✔
1251
  try {
1252
    model::cells[index]->set_temperature(T, instance_index, set_contained);
80!
1253
  } catch (const std::exception& e) {
×
1254
    set_errmsg(e.what());
×
1255
    return OPENMC_E_UNASSIGNED;
×
1256
  }
×
1257
  return 0;
80✔
1258
}
1259

1260
extern "C" int openmc_cell_set_density(
80✔
1261
  int32_t index, double density, const int32_t* instance, bool set_contained)
1262
{
1263
  if (index < 0 || index >= model::cells.size()) {
80!
1264
    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
×
1265
    return OPENMC_E_OUT_OF_BOUNDS;
×
1266
  }
1267

1268
  int32_t instance_index = instance ? *instance : -1;
80✔
1269
  try {
1270
    model::cells[index]->set_density(density, instance_index, set_contained);
80!
1271
  } catch (const std::exception& e) {
×
1272
    set_errmsg(e.what());
×
1273
    return OPENMC_E_UNASSIGNED;
×
1274
  }
×
1275
  return 0;
80✔
1276
}
1277

1278
extern "C" int openmc_cell_get_temperature(
80✔
1279
  int32_t index, const int32_t* instance, double* T)
1280
{
1281
  if (index < 0 || index >= model::cells.size()) {
80!
1282
    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
×
1283
    return OPENMC_E_OUT_OF_BOUNDS;
×
1284
  }
1285

1286
  int32_t instance_index = instance ? *instance : -1;
80✔
1287
  try {
1288
    *T = model::cells[index]->temperature(instance_index);
80!
1289
  } catch (const std::exception& e) {
×
1290
    set_errmsg(e.what());
×
1291
    return OPENMC_E_UNASSIGNED;
×
1292
  }
×
1293
  return 0;
80✔
1294
}
1295

1296
extern "C" int openmc_cell_get_density(
80✔
1297
  int32_t index, const int32_t* instance, double* density)
1298
{
1299
  if (index < 0 || index >= model::cells.size()) {
80!
1300
    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
×
1301
    return OPENMC_E_OUT_OF_BOUNDS;
×
1302
  }
1303

1304
  int32_t instance_index = instance ? *instance : -1;
80✔
1305
  try {
1306
    if (model::cells[index]->type_ != Fill::MATERIAL) {
80!
1307
      fatal_error(
×
1308
        fmt::format("Cell {}, instance {} is not filled with a material.",
×
1309
          model::cells[index]->id_, instance_index));
×
1310
    }
1311

1312
    int32_t mat_index = model::cells[index]->material(instance_index);
80✔
1313
    if (mat_index == MATERIAL_VOID) {
80!
1314
      *density = 0.0;
×
1315
    } else {
1316
      *density = model::cells[index]->density_mult(instance_index) *
160!
1317
                 model::materials[mat_index]->density_gpcc();
80✔
1318
    }
1319
  } catch (const std::exception& e) {
×
1320
    set_errmsg(e.what());
×
1321
    return OPENMC_E_UNASSIGNED;
×
1322
  }
×
1323
  return 0;
80✔
1324
}
1325

1326
//! Get the bounding box of a cell
1327
extern "C" int openmc_cell_bounding_box(
50✔
1328
  const int32_t index, double* llc, double* urc)
1329
{
1330

1331
  BoundingBox bbox;
50✔
1332

1333
  const auto& c = model::cells[index];
50✔
1334
  bbox = c->bounding_box();
50!
1335

1336
  // set lower left corner values
1337
  llc[0] = bbox.xmin;
50✔
1338
  llc[1] = bbox.ymin;
50✔
1339
  llc[2] = bbox.zmin;
50✔
1340

1341
  // set upper right corner values
1342
  urc[0] = bbox.xmax;
50✔
1343
  urc[1] = bbox.ymax;
50✔
1344
  urc[2] = bbox.zmax;
50✔
1345

1346
  return 0;
50✔
1347
}
1348

1349
//! Get the name of a cell
1350
extern "C" int openmc_cell_get_name(int32_t index, const char** name)
20✔
1351
{
1352
  if (index < 0 || index >= model::cells.size()) {
20!
1353
    set_errmsg("Index in cells array is out of bounds.");
×
1354
    return OPENMC_E_OUT_OF_BOUNDS;
×
1355
  }
1356

1357
  *name = model::cells[index]->name().data();
20✔
1358

1359
  return 0;
20✔
1360
}
1361

1362
//! Set the name of a cell
1363
extern "C" int openmc_cell_set_name(int32_t index, const char* name)
10✔
1364
{
1365
  if (index < 0 || index >= model::cells.size()) {
10!
1366
    set_errmsg("Index in cells array is out of bounds.");
×
1367
    return OPENMC_E_OUT_OF_BOUNDS;
×
1368
  }
1369

1370
  model::cells[index]->set_name(name);
10!
1371

1372
  return 0;
10✔
1373
}
1374

1375
//==============================================================================
1376
//! Define a containing (parent) cell
1377
//==============================================================================
1378

1379
//! Used to locate a universe fill in the geometry
1380
struct ParentCell {
1381
  bool operator==(const ParentCell& other) const
126✔
1382
  {
1383
    return cell_index == other.cell_index &&
252!
1384
           lattice_index == other.lattice_index;
252!
1385
  }
1386

1387
  bool operator<(const ParentCell& other) const
1388
  {
1389
    return cell_index < other.cell_index ||
1390
           (cell_index == other.cell_index &&
1391
             lattice_index < other.lattice_index);
1392
  }
1393

1394
  int64_t cell_index;
1395
  int64_t lattice_index;
1396
};
1397

1398
//! Structure used to insert ParentCell into hashed STL data structures
1399
struct ParentCellHash {
1400
  std::size_t operator()(const ParentCell& p) const
588✔
1401
  {
1402
    return 4096 * p.cell_index + p.lattice_index;
588✔
1403
  }
1404
};
1405

1406
//! Used to manage a traversal stack when locating parent cells of a cell
1407
//! instance in the model
1408
struct ParentCellStack {
1409

1410
  //! push method that adds to the parent_cells visited cells for this search
1411
  //! universe
1412
  void push(int32_t search_universe, const ParentCell& pc)
98✔
1413
  {
1414
    parent_cells_.push_back(pc);
98✔
1415
    // add parent cell to the set of cells we've visited for this search
1416
    // universe
1417
    visited_cells_[search_universe].insert(pc);
98✔
1418
  }
98✔
1419

1420
  //! removes the last parent_cell and clears the visited cells for the popped
1421
  //! cell's universe
1422
  void pop()
70✔
1423
  {
1424
    visited_cells_[this->current_univ()].clear();
70!
1425
    parent_cells_.pop_back();
70✔
1426
  }
70✔
1427

1428
  //! checks whether or not the parent cell has been visited already for this
1429
  //! search universe
1430
  bool visited(int32_t search_universe, const ParentCell& parent_cell)
490✔
1431
  {
1432
    return visited_cells_[search_universe].count(parent_cell) != 0;
490✔
1433
  }
1434

1435
  //! return the next universe to search for a parent cell
1436
  int32_t current_univ() const
70✔
1437
  {
1438
    return model::cells[parent_cells_.back().cell_index]->universe_;
70✔
1439
  }
1440

1441
  //! indicates whether nor not parent cells are present on the stack
1442
  bool empty() const { return parent_cells_.empty(); }
70✔
1443

1444
  //! compute an instance for the provided distribcell index
1445
  int32_t compute_instance(int32_t distribcell_index) const
168✔
1446
  {
1447
    if (distribcell_index == C_NONE)
168✔
1448
      return 0;
56✔
1449

1450
    int32_t instance = 0;
112✔
1451
    for (const auto& parent_cell : this->parent_cells_) {
210✔
1452
      auto& cell = model::cells[parent_cell.cell_index];
98✔
1453
      if (cell->type_ == Fill::UNIVERSE) {
98!
1454
        instance += cell->offset_[distribcell_index];
×
1455
      } else if (cell->type_ == Fill::LATTICE) {
98!
1456
        auto& lattice = model::lattices[cell->fill_];
98✔
1457
        instance +=
98✔
1458
          lattice->offset(distribcell_index, parent_cell.lattice_index);
98!
1459
      }
1460
    }
1461
    return instance;
112✔
1462
  }
1463

1464
  // Accessors
1465
  vector<ParentCell>& parent_cells() { return parent_cells_; }
294✔
1466
  const vector<ParentCell>& parent_cells() const { return parent_cells_; }
1467

1468
  // Data Members
1469
  vector<ParentCell> parent_cells_;
1470
  std::unordered_map<int32_t, std::unordered_set<ParentCell, ParentCellHash>>
1471
    visited_cells_;
1472
};
1473

1474
vector<ParentCell> Cell::find_parent_cells(
×
1475
  int32_t instance, const Position& r) const
1476
{
1477

1478
  // create a temporary particle
1479
  GeometryState dummy_particle {};
×
1480
  dummy_particle.r() = r;
×
1481
  dummy_particle.u() = {0., 0., 1.};
×
1482

1483
  return find_parent_cells(instance, dummy_particle);
×
1484
}
×
1485

1486
vector<ParentCell> Cell::find_parent_cells(
×
1487
  int32_t instance, GeometryState& p) const
1488
{
1489
  // look up the particle's location
1490
  exhaustive_find_cell(p);
×
1491
  const auto& coords = p.coord();
×
1492

1493
  // build a parent cell stack from the particle coordinates
1494
  ParentCellStack stack;
×
1495
  bool cell_found = false;
×
1496
  for (auto it = coords.begin(); it != coords.end(); it++) {
×
1497
    const auto& coord = *it;
×
1498
    const auto& cell = model::cells[coord.cell()];
×
1499
    // if the cell at this level matches the current cell, stop adding to the
1500
    // stack
1501
    if (coord.cell() == model::cell_map[this->id_]) {
×
1502
      cell_found = true;
×
1503
      break;
×
1504
    }
1505

1506
    // if filled with a lattice, get the lattice index from the next
1507
    // level in the coordinates to push to the stack
1508
    int lattice_idx = C_NONE;
×
1509
    if (cell->type_ == Fill::LATTICE) {
×
1510
      const auto& next_coord = *(it + 1);
×
1511
      lattice_idx = model::lattices[next_coord.lattice()]->get_flat_index(
×
1512
        next_coord.lattice_index());
1513
    }
1514
    stack.push(coord.universe(), {coord.cell(), lattice_idx});
×
1515
  }
1516

1517
  // if this loop finished because the cell was found and
1518
  // the instance matches the one requested in the call
1519
  // we have the correct path and can return the stack
1520
  if (cell_found &&
×
1521
      stack.compute_instance(this->distribcell_index_) == instance) {
×
1522
    return stack.parent_cells();
×
1523
  }
1524

1525
  // fall back on an exhaustive search for the cell's parents
1526
  return exhaustive_find_parent_cells(instance);
×
1527
}
×
1528

1529
vector<ParentCell> Cell::exhaustive_find_parent_cells(int32_t instance) const
98✔
1530
{
1531
  ParentCellStack stack;
98✔
1532
  // start with this cell's universe
1533
  int32_t prev_univ_idx;
1534
  int32_t univ_idx = this->universe_;
98✔
1535

1536
  while (true) {
1537
    const auto& univ = model::universes[univ_idx];
168✔
1538
    prev_univ_idx = univ_idx;
168✔
1539

1540
    // search for a cell that is filled w/ this universe
1541
    for (const auto& cell : model::cells) {
1,078✔
1542
      // if this is a material-filled cell, move on
1543
      if (cell->type_ == Fill::MATERIAL)
1,008✔
1544
        continue;
560✔
1545

1546
      if (cell->type_ == Fill::UNIVERSE) {
448✔
1547
        // if this is in the set of cells previously visited for this universe,
1548
        // move on
1549
        if (stack.visited(univ_idx, {model::cell_map[cell->id_], C_NONE}))
266!
1550
          continue;
×
1551

1552
        // if this cell contains the universe we're searching for, add it to the
1553
        // stack
1554
        if (cell->fill_ == univ_idx) {
266!
1555
          stack.push(univ_idx, {model::cell_map[cell->id_], C_NONE});
×
1556
          univ_idx = cell->universe_;
×
1557
        }
1558
      } else if (cell->type_ == Fill::LATTICE) {
182!
1559
        // retrieve the lattice and lattice universes
1560
        const auto& lattice = model::lattices[cell->fill_];
182✔
1561
        const auto& lattice_univs = lattice->universes_;
182✔
1562

1563
        // start search for universe
1564
        auto lat_it = lattice_univs.begin();
182✔
1565
        while (true) {
1566
          // find the next lattice cell with this universe
1567
          lat_it = std::find(lat_it, lattice_univs.end(), univ_idx);
308!
1568
          if (lat_it == lattice_univs.end())
308✔
1569
            break;
84✔
1570

1571
          int lattice_idx = lat_it - lattice_univs.begin();
224✔
1572

1573
          // move iterator forward one to avoid finding the same entry
1574
          lat_it++;
224✔
1575
          if (stack.visited(
224✔
1576
                univ_idx, {model::cell_map[cell->id_], lattice_idx}))
224!
1577
            continue;
126✔
1578

1579
          // add this cell and lattice index to the stack and exit loop
1580
          stack.push(univ_idx, {model::cell_map[cell->id_], lattice_idx});
98!
1581
          univ_idx = cell->universe_;
98✔
1582
          break;
98✔
1583
        }
126✔
1584
      }
1585
      // if we've updated the universe, break
1586
      if (prev_univ_idx != univ_idx)
448✔
1587
        break;
98✔
1588
    } // end cell loop search for universe
1589

1590
    // if we're at the top of the geometry and the instance matches, we're done
1591
    if (univ_idx == model::root_universe &&
336!
1592
        stack.compute_instance(this->distribcell_index_) == instance)
168!
1593
      break;
98✔
1594

1595
    // if there is no match on the original cell's universe, report an error
1596
    if (univ_idx == this->universe_) {
70!
1597
      fatal_error(
×
1598
        fmt::format("Could not find the parent cells for cell {}, instance {}.",
×
1599
          this->id_, instance));
×
1600
    }
1601

1602
    // if we don't find a suitable update, adjust the stack and continue
1603
    if (univ_idx == model::root_universe || univ_idx == prev_univ_idx) {
70!
1604
      stack.pop();
70!
1605
      univ_idx = stack.empty() ? this->universe_ : stack.current_univ();
70!
1606
    }
1607

1608
  } // end while
70✔
1609

1610
  // reverse the stack so the highest cell comes first
1611
  std::reverse(stack.parent_cells().begin(), stack.parent_cells().end());
98!
1612
  return stack.parent_cells();
196!
1613
}
98✔
1614

1615
std::unordered_map<int32_t, vector<int32_t>> Cell::get_contained_cells(
140✔
1616
  int32_t instance, Position* hint) const
1617
{
1618
  std::unordered_map<int32_t, vector<int32_t>> contained_cells;
140✔
1619

1620
  // if this is a material-filled cell it has no contained cells
1621
  if (this->type_ == Fill::MATERIAL)
140✔
1622
    return contained_cells;
42✔
1623

1624
  // find the pathway through the geometry to this cell
1625
  vector<ParentCell> parent_cells;
98✔
1626

1627
  // if a positional hint is provided, attempt to do a fast lookup
1628
  // of the parent cells
1629
  parent_cells = hint ? find_parent_cells(instance, *hint)
196!
1630
                      : exhaustive_find_parent_cells(instance);
98✔
1631

1632
  // if this cell is filled w/ a material, it contains no other cells
1633
  if (type_ != Fill::MATERIAL) {
98!
1634
    this->get_contained_cells_inner(contained_cells, parent_cells);
98!
1635
  }
1636

1637
  return contained_cells;
98✔
1638
}
98✔
1639

1640
//! Get all cells within this cell
1641
void Cell::get_contained_cells_inner(
76,790✔
1642
  std::unordered_map<int32_t, vector<int32_t>>& contained_cells,
1643
  vector<ParentCell>& parent_cells) const
1644
{
1645

1646
  // filled by material, determine instance based on parent cells
1647
  if (type_ == Fill::MATERIAL) {
76,790✔
1648
    int instance = 0;
76,244✔
1649
    if (this->distribcell_index_ >= 0) {
76,244!
1650
      for (auto& parent_cell : parent_cells) {
228,732✔
1651
        auto& cell = model::cells[parent_cell.cell_index];
152,488✔
1652
        if (cell->type_ == Fill::UNIVERSE) {
152,488✔
1653
          instance += cell->offset_[distribcell_index_];
74,844✔
1654
        } else if (cell->type_ == Fill::LATTICE) {
77,644!
1655
          auto& lattice = model::lattices[cell->fill_];
77,644✔
1656
          instance += lattice->offset(
155,288✔
1657
            this->distribcell_index_, parent_cell.lattice_index);
77,644!
1658
        }
1659
      }
1660
    }
1661
    // add entry to contained cells
1662
    contained_cells[model::cell_map[id_]].push_back(instance);
76,244!
1663
    // filled with universe, add the containing cell to the parent cells
1664
    // and recurse
1665
  } else if (type_ == Fill::UNIVERSE) {
546✔
1666
    parent_cells.push_back({model::cell_map[id_], -1});
462!
1667
    auto& univ = model::universes[fill_];
462✔
1668
    for (auto cell_index : univ->cells_) {
2,772✔
1669
      auto& cell = model::cells[cell_index];
2,310✔
1670
      cell->get_contained_cells_inner(contained_cells, parent_cells);
2,310!
1671
    }
1672
    parent_cells.pop_back();
462✔
1673
    // filled with a lattice, visit each universe in the lattice
1674
    // with a recursive call to collect the cell instances
1675
  } else if (type_ == Fill::LATTICE) {
84!
1676
    auto& lattice = model::lattices[fill_];
84✔
1677
    for (auto i = lattice->begin(); i != lattice->end(); ++i) {
74,172!
1678
      auto& univ = model::universes[*i];
74,088✔
1679
      parent_cells.push_back({model::cell_map[id_], i.indx_});
74,088!
1680
      for (auto cell_index : univ->cells_) {
148,470✔
1681
        auto& cell = model::cells[cell_index];
74,382✔
1682
        cell->get_contained_cells_inner(contained_cells, parent_cells);
74,382!
1683
      }
1684
      parent_cells.pop_back();
74,088✔
1685
    }
1686
  }
1687
}
76,790✔
1688

1689
//! Return the index in the cells array of a cell with a given ID
1690
extern "C" int openmc_get_cell_index(int32_t id, int32_t* index)
400✔
1691
{
1692
  auto it = model::cell_map.find(id);
400!
1693
  if (it != model::cell_map.end()) {
400✔
1694
    *index = it->second;
390✔
1695
    return 0;
390✔
1696
  } else {
1697
    set_errmsg("No cell exists with ID=" + std::to_string(id) + ".");
10!
1698
    return OPENMC_E_INVALID_ID;
10✔
1699
  }
1700
}
1701

1702
//! Return the ID of a cell
1703
extern "C" int openmc_cell_get_id(int32_t index, int32_t* id)
400,840✔
1704
{
1705
  if (index >= 0 && index < model::cells.size()) {
400,840!
1706
    *id = model::cells[index]->id_;
400,840✔
1707
    return 0;
400,840✔
1708
  } else {
1709
    set_errmsg("Index in cells array is out of bounds.");
×
1710
    return OPENMC_E_OUT_OF_BOUNDS;
×
1711
  }
1712
}
1713

1714
//! Set the ID of a cell
1715
extern "C" int openmc_cell_set_id(int32_t index, int32_t id)
20✔
1716
{
1717
  if (index >= 0 && index < model::cells.size()) {
20!
1718
    model::cells[index]->id_ = id;
20✔
1719
    model::cell_map[id] = index;
20✔
1720
    return 0;
20✔
1721
  } else {
1722
    set_errmsg("Index in cells array is out of bounds.");
×
1723
    return OPENMC_E_OUT_OF_BOUNDS;
×
1724
  }
1725
}
1726

1727
//! Return the translation vector of a cell
1728
extern "C" int openmc_cell_get_translation(int32_t index, double xyz[])
50✔
1729
{
1730
  if (index >= 0 && index < model::cells.size()) {
50!
1731
    auto& cell = model::cells[index];
50✔
1732
    xyz[0] = cell->translation_.x;
50✔
1733
    xyz[1] = cell->translation_.y;
50✔
1734
    xyz[2] = cell->translation_.z;
50✔
1735
    return 0;
50✔
1736
  } else {
1737
    set_errmsg("Index in cells array is out of bounds.");
×
1738
    return OPENMC_E_OUT_OF_BOUNDS;
×
1739
  }
1740
}
1741

1742
//! Set the translation vector of a cell
1743
extern "C" int openmc_cell_set_translation(int32_t index, const double xyz[])
30✔
1744
{
1745
  if (index >= 0 && index < model::cells.size()) {
30!
1746
    if (model::cells[index]->fill_ == C_NONE) {
30✔
1747
      set_errmsg(fmt::format("Cannot apply a translation to cell {}"
10!
1748
                             " because it is not filled with another universe",
1749
        index));
1750
      return OPENMC_E_GEOMETRY;
10✔
1751
    }
1752
    model::cells[index]->translation_ = Position(xyz);
20✔
1753
    return 0;
20✔
1754
  } else {
1755
    set_errmsg("Index in cells array is out of bounds.");
×
1756
    return OPENMC_E_OUT_OF_BOUNDS;
×
1757
  }
1758
}
1759

1760
//! Return the rotation matrix of a cell
1761
extern "C" int openmc_cell_get_rotation(int32_t index, double rot[], size_t* n)
50✔
1762
{
1763
  if (index >= 0 && index < model::cells.size()) {
50!
1764
    auto& cell = model::cells[index];
50✔
1765
    *n = cell->rotation_.size();
50✔
1766
    std::memcpy(rot, cell->rotation_.data(), *n * sizeof(cell->rotation_[0]));
50✔
1767
    return 0;
50✔
1768
  } else {
1769
    set_errmsg("Index in cells array is out of bounds.");
×
1770
    return OPENMC_E_OUT_OF_BOUNDS;
×
1771
  }
1772
}
1773

1774
//! Set the flattened rotation matrix of a cell
1775
extern "C" int openmc_cell_set_rotation(
40✔
1776
  int32_t index, const double rot[], size_t rot_len)
1777
{
1778
  if (index >= 0 && index < model::cells.size()) {
40!
1779
    if (model::cells[index]->fill_ == C_NONE) {
40✔
1780
      set_errmsg(fmt::format("Cannot apply a rotation to cell {}"
10!
1781
                             " because it is not filled with another universe",
1782
        index));
1783
      return OPENMC_E_GEOMETRY;
10✔
1784
    }
1785
    std::vector<double> vec_rot(rot, rot + rot_len);
30!
1786
    model::cells[index]->set_rotation(vec_rot);
30!
1787
    return 0;
30✔
1788
  } else {
30✔
1789
    set_errmsg("Index in cells array is out of bounds.");
×
1790
    return OPENMC_E_OUT_OF_BOUNDS;
×
1791
  }
1792
}
1793

1794
//! Get the number of instances of the requested cell
1795
extern "C" int openmc_cell_get_num_instances(
10✔
1796
  int32_t index, int32_t* num_instances)
1797
{
1798
  if (index < 0 || index >= model::cells.size()) {
10!
1799
    set_errmsg("Index in cells array is out of bounds.");
×
1800
    return OPENMC_E_OUT_OF_BOUNDS;
×
1801
  }
1802
  *num_instances = model::cells[index]->n_instances();
10✔
1803
  return 0;
10✔
1804
}
1805

1806
//! Extend the cells array by n elements
1807
extern "C" int openmc_extend_cells(
20✔
1808
  int32_t n, int32_t* index_start, int32_t* index_end)
1809
{
1810
  if (index_start)
20!
1811
    *index_start = model::cells.size();
20✔
1812
  if (index_end)
20!
1813
    *index_end = model::cells.size() + n - 1;
×
1814
  for (int32_t i = 0; i < n; i++) {
40✔
1815
    model::cells.push_back(make_unique<CSGCell>());
20!
1816
  }
1817
  return 0;
20✔
1818
}
1819

1820
extern "C" int cells_size()
40✔
1821
{
1822
  return model::cells.size();
40✔
1823
}
1824

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

© 2025 Coveralls, Inc