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

openmc-dev / openmc / 14840340664

05 May 2025 03:38PM UTC coverage: 85.195% (-0.009%) from 85.204%
14840340664

Pull #3392

github

web-flow
Merge 5bc1ec35f into 1e7d8324e
Pull Request #3392: Map Compton subshell data to atomic relaxation data

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

330 existing lines in 19 files now uncovered.

52194 of 61264 relevant lines covered (85.2%)

37398320.1 hits per line

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

75.0
/include/openmc/cell.h
1
#ifndef OPENMC_CELL_H
2
#define OPENMC_CELL_H
3

4
#include <cstdint>
5
#include <functional> // for hash
6
#include <limits>
7
#include <string>
8
#include <unordered_map>
9
#include <unordered_set>
10

11
#include "hdf5.h"
12
#include "pugixml.hpp"
13

14
#include "openmc/bounding_box.h"
15
#include "openmc/constants.h"
16
#include "openmc/memory.h" // for unique_ptr
17
#include "openmc/neighbor_list.h"
18
#include "openmc/position.h"
19
#include "openmc/surface.h"
20
#include "openmc/universe.h"
21
#include "openmc/vector.h"
22

23
namespace openmc {
24

25
//==============================================================================
26
// Constants
27
//==============================================================================
28

29
enum class Fill { MATERIAL, UNIVERSE, LATTICE };
30

31
constexpr int32_t OP_LEFT_PAREN {std::numeric_limits<int32_t>::max()};
32
constexpr int32_t OP_RIGHT_PAREN {std::numeric_limits<int32_t>::max() - 1};
33
constexpr int32_t OP_COMPLEMENT {std::numeric_limits<int32_t>::max() - 2};
34
constexpr int32_t OP_INTERSECTION {std::numeric_limits<int32_t>::max() - 3};
35
constexpr int32_t OP_UNION {std::numeric_limits<int32_t>::max() - 4};
36

37
//==============================================================================
38
// Global variables
39
//==============================================================================
40

41
class Cell;
42
class GeometryState;
43
class ParentCell;
44
class CellInstance;
45
class Universe;
46
class UniversePartitioner;
47

48
namespace model {
49
extern std::unordered_map<int32_t, int32_t> cell_map;
50
extern vector<unique_ptr<Cell>> cells;
51

52
} // namespace model
53

54
//==============================================================================
55

56
class Region {
57
public:
58
  //----------------------------------------------------------------------------
59
  // Constructors
60
  Region() {}
30,904✔
61
  explicit Region(std::string region_spec, int32_t cell_id);
62

63
  //----------------------------------------------------------------------------
64
  // Methods
65

66
  //! \brief Determine if a cell contains the particle at a given location.
67
  //!
68
  //! The bounds of the cell are determined by a logical expression involving
69
  //! surface half-spaces. The expression used is given in infix notation
70
  //!
71
  //! The function is split into two cases, one for simple cells (those
72
  //! involving only the intersection of half-spaces) and one for complex cells.
73
  //! Both cases use short circuiting; however, in the case fo complex cells,
74
  //! the complexity increases with the binary operators involved.
75
  //! \param r The 3D Cartesian coordinate to check.
76
  //! \param u A direction used to "break ties" the coordinates are very
77
  //!   close to a surface.
78
  //! \param on_surface The signed index of a surface that the coordinate is
79
  //!   known to be on.  This index takes precedence over surface sense
80
  //!   calculations.
81
  bool contains(Position r, Direction u, int32_t on_surface) const;
82

83
  //! Find the oncoming boundary of this cell.
84
  std::pair<double, int32_t> distance(
85
    Position r, Direction u, int32_t on_surface) const;
86

87
  //! Get the BoundingBox for this cell.
88
  BoundingBox bounding_box(int32_t cell_id) const;
89

90
  //! Get the CSG expression as a string
91
  std::string str() const;
92

93
  //! Get a vector containing all the surfaces in the region expression
94
  vector<int32_t> surfaces() const;
95

96
  //----------------------------------------------------------------------------
97
  // Accessors
98

99
  //! Get Boolean of if the cell is simple or not
100
  bool is_simple() const { return simple_; }
2,147,483,647✔
101

102
private:
103
  //----------------------------------------------------------------------------
104
  // Private Methods
105

106
  //! Get a vector of the region expression in postfix notation
107
  vector<int32_t> generate_postfix(int32_t cell_id) const;
108

109
  //! Determine if a particle is inside the cell for a simple cell (only
110
  //! intersection operators)
111
  bool contains_simple(Position r, Direction u, int32_t on_surface) const;
112

113
  //! Determine if a particle is inside the cell for a complex cell.
114
  //!
115
  //! Uses the comobination of half-spaces and binary operators to determine
116
  //! if short circuiting can be used. Short cicuiting uses the relative and
117
  //! absolute depth of parentheses in the expression.
118
  bool contains_complex(Position r, Direction u, int32_t on_surface) const;
119

120
  //! BoundingBox if the paritcle is in a simple cell.
121
  BoundingBox bounding_box_simple() const;
122

123
  //! BoundingBox if the particle is in a complex cell.
124
  BoundingBox bounding_box_complex(vector<int32_t> postfix) const;
125

126
  //! Enfource precedence: Parenthases, Complement, Intersection, Union
127
  void add_precedence();
128

129
  //! Add parenthesis to enforce precedence
130
  int64_t add_parentheses(int64_t start);
131

132
  //! Remove complement operators from the expression
133
  void remove_complement_ops();
134

135
  //! Remove complement operators by using DeMorgan's laws
136
  void apply_demorgan(
137
    vector<int32_t>::iterator start, vector<int32_t>::iterator stop);
138

139
  //----------------------------------------------------------------------------
140
  // Private Data
141

142
  //! Definition of spatial region as Boolean expression of half-spaces
143
  // TODO: Should this be a vector of some other type
144
  vector<int32_t> expression_;
145
  bool simple_; //!< Does the region contain only intersections?
146
};
147

148
//==============================================================================
149

150
class Cell {
151
public:
152
  //----------------------------------------------------------------------------
153
  // Constructors, destructors, factory functions
154

155
  explicit Cell(pugi::xml_node cell_node);
156
  Cell() {};
31,100✔
157
  virtual ~Cell() = default;
31,098✔
UNCOV
158

×
159
  //----------------------------------------------------------------------------
31,098✔
160
  // Methods
161

162
  //! \brief Determine if a cell contains the particle at a given location.
163
  //!
164
  //! The bounds of the cell are detemined by a logical expression involving
165
  //! surface half-spaces. At initialization, the expression was converted
166
  //! to RPN notation.
167
  //!
168
  //! The function is split into two cases, one for simple cells (those
169
  //! involving only the intersection of half-spaces) and one for complex cells.
170
  //! Simple cells can be evaluated with short circuit evaluation, i.e., as soon
171
  //! as we know that one half-space is not satisfied, we can exit. This
172
  //! provides a performance benefit for the common case. In
173
  //! contains_complex, we evaluate the RPN expression using a stack, similar to
174
  //! how a RPN calculator would work.
175
  //! \param r The 3D Cartesian coordinate to check.
176
  //! \param u A direction used to "break ties" the coordinates are very
177
  //!   close to a surface.
178
  //! \param on_surface The signed index of a surface that the coordinate is
179
  //!   known to be on.  This index takes precedence over surface sense
180
  //!   calculations.
181
  virtual bool contains(Position r, Direction u, int32_t on_surface) const = 0;
182

183
  //! Find the oncoming boundary of this cell.
184
  virtual std::pair<double, int32_t> distance(
185
    Position r, Direction u, int32_t on_surface, GeometryState* p) const = 0;
186

187
  //! Write all information needed to reconstruct the cell to an HDF5 group.
188
  //! \param group_id An HDF5 group id.
189
  void to_hdf5(hid_t group_id) const;
190

191
  virtual void to_hdf5_inner(hid_t group_id) const = 0;
192

193
  //! Export physical properties to HDF5
194
  //! \param[in] group  HDF5 group to read from
195
  void export_properties_hdf5(hid_t group) const;
196

197
  //! Import physical properties from HDF5
198
  //! \param[in] group  HDF5 group to write to
199
  void import_properties_hdf5(hid_t group);
200

201
  //! Get the BoundingBox for this cell.
202
  virtual BoundingBox bounding_box() const = 0;
203

204
  //! Get a vector of surfaces in the cell
205
  virtual vector<int32_t> surfaces() const { return vector<int32_t>(); }
206

UNCOV
207
  //! Check if the cell region expression is simple
×
208
  virtual bool is_simple() const { return true; }
209

210
  //----------------------------------------------------------------------------
92,429✔
211
  // Accessors
212

213
  //! Get the temperature of a cell instance
214
  //! \param[in] instance Instance index. If -1 is given, the temperature for
215
  //!   the first instance is returned.
216
  //! \return Temperature in [K]
217
  double temperature(int32_t instance = -1) const;
218

219
  //! Set the temperature of a cell instance
220
  //! \param[in] T Temperature in [K]
221
  //! \param[in] instance Instance index. If -1 is given, the temperature for
222
  //!   all instances is set.
223
  //! \param[in] set_contained If this cell is not filled with a material,
224
  //!   collect all contained cells with material fills and set their
225
  //!   temperatures.
226
  void set_temperature(
227
    double T, int32_t instance = -1, bool set_contained = false);
228

229
  //! Set the rotation matrix of a cell instance
230
  //! \param[in] rot The rotation matrix of length 3 or 9
231
  void set_rotation(const vector<double>& rot);
232

233
  //! Get the name of a cell
234
  //! \return Cell name
235
  const std::string& name() const { return name_; };
236

237
  //! Set the temperature of a cell instance
288✔
238
  //! \param[in] name Cell name
239
  void set_name(const std::string& name) { name_ = name; };
240

241
  //! Get all cell instances contained by this cell
11✔
242
  //! \param[in] instance Instance of the cell for which to get contained cells
243
  //! (default instance is zero)
244
  //! \param[in] hint positional hint for determining the parent cells
245
  //! \return Map with cell indexes as keys and
246
  //! instances as values
247
  std::unordered_map<int32_t, vector<int32_t>> get_contained_cells(
248
    int32_t instance = 0, Position* hint = nullptr) const;
249

250
  //! Determine the material index corresponding to a specific cell instance,
251
  //! taking into account presence of distribcell material
252
  //! \param[in] instance of the cell
253
  //! \return material index
254
  int32_t material(int32_t instance) const
×
255
  {
256
    // If distributed materials are used, then each instance has its own
257
    // material definition. If distributed materials are not used, then
258
    // all instances used the same material stored at material_[0]. The
259
    // presence of distributed materials is inferred from the size of
260
    // the material_ vector being greater than one.
261
    if (material_.size() > 1) {
×
262
      return material_[instance];
×
263
    } else {
264
      return material_[0];
×
265
    }
266
  }
267

268
  //! Determine the temperature index corresponding to a specific cell instance,
269
  //! taking into account presence of distribcell temperature
270
  //! \param[in] instance of the cell
271
  //! \return temperature index
272
  double sqrtkT(int32_t instance) const
273
  {
274
    // If distributed materials are used, then each instance has its own
275
    // temperature definition. If distributed materials are not used, then
276
    // all instances used the same temperature stored at sqrtkT_[0]. The
277
    // presence of distributed materials is inferred from the size of
278
    // the sqrtkT_ vector being greater than one.
279
    if (sqrtkT_.size() > 1) {
280
      return sqrtkT_[instance];
281
    } else {
282
      return sqrtkT_[0];
283
    }
284
  }
285

286
protected:
287
  //! Determine the path to this cell instance in the geometry hierarchy
288
  //! \param[in] instance of the cell to find parent cells for
289
  //! \param[in] r position used to do a fast search for parent cells
290
  //! \return parent cells
291
  vector<ParentCell> find_parent_cells(
292
    int32_t instance, const Position& r) const;
293

294
  //! Determine the path to this cell instance in the geometry hierarchy
295
  //! \param[in] instance of the cell to find parent cells for
296
  //! \param[in] p particle used to do a fast search for parent cells
297
  //! \return parent cells
298
  vector<ParentCell> find_parent_cells(
299
    int32_t instance, GeometryState& p) const;
300

301
  //! Determine the path to this cell instance in the geometry hierarchy
302
  //! \param[in] instance of the cell to find parent cells for
303
  //! \return parent cells
304
  vector<ParentCell> exhaustive_find_parent_cells(int32_t instance) const;
305

306
  //! Inner function for retrieving contained cells
307
  void get_contained_cells_inner(
308
    std::unordered_map<int32_t, vector<int32_t>>& contained_cells,
309
    vector<ParentCell>& parent_cells) const;
310

311
public:
312
  //----------------------------------------------------------------------------
313
  // Data members
314

315
  int32_t id_;              //!< Unique ID
316
  std::string name_;        //!< User-defined name
317
  Fill type_;               //!< Material, universe, or lattice
318
  int32_t universe_;        //!< Universe # this cell is in
319
  int32_t fill_;            //!< Universe # filling this cell
320
  int32_t n_instances_ {0}; //!< Number of instances of this cell
321

322
  //! \brief Index corresponding to this cell in distribcell arrays
323
  int distribcell_index_ {C_NONE};
324

325
  //! \brief Material(s) within this cell.
326
  //!
327
  //! May be multiple materials for distribcell.
328
  vector<int32_t> material_;
329

330
  //! \brief Temperature(s) within this cell.
331
  //!
332
  //! The stored values are actually sqrt(k_Boltzmann * T) for each temperature
333
  //! T. The units are sqrt(eV).
334
  vector<double> sqrtkT_;
335

336
  //! \brief Neighboring cells in the same universe.
337
  NeighborList neighbors_;
338

339
  Position translation_ {0, 0, 0}; //!< Translation vector for filled universe
340

341
  //! \brief Rotational tranfsormation of the filled universe.
342
  //
343
  //! The vector is empty if there is no rotation. Otherwise, the first 9 values
344
  //! give the rotation matrix in row-major order. When the user specifies
345
  //! rotation angles about the x-, y- and z- axes in degrees, these values are
346
  //! also present at the end of the vector, making it of length 12.
347
  vector<double> rotation_;
348

349
  vector<int32_t> offset_; //!< Distribcell offset table
350

351
  // Right now, either CSG or DAGMC cells are used.
352
  virtual GeometryType geom_type() const = 0;
353
};
354

355
struct CellInstanceItem {
356
  int32_t index {-1};    //! Index into global cells array
357
  int lattice_indx {-1}; //! Flat index value of the lattice cell
358
};
359

360
//==============================================================================
361

362
class CSGCell : public Cell {
363
public:
364
  //----------------------------------------------------------------------------
365
  // Constructors
366
  CSGCell() = default;
367
  explicit CSGCell(pugi::xml_node cell_node);
368

22✔
369
  //----------------------------------------------------------------------------
370
  // Methods
371
  vector<int32_t> surfaces() const override { return region_.surfaces(); }
372

373
  std::pair<double, int32_t> distance(Position r, Direction u,
5,389✔
374
    int32_t on_surface, GeometryState* p) const override
375
  {
2,147,483,647✔
376
    return region_.distance(r, u, on_surface);
377
  }
378

2,147,483,647✔
379
  bool contains(Position r, Direction u, int32_t on_surface) const override
380
  {
381
    return region_.contains(r, u, on_surface);
2,147,483,647✔
382
  }
383

2,147,483,647✔
384
  BoundingBox bounding_box() const override
385
  {
386
    return region_.bounding_box(id_);
88✔
387
  }
388

88✔
389
  void to_hdf5_inner(hid_t group_id) const override;
390

391
  bool is_simple() const override { return region_.is_simple(); }
392

393
  virtual GeometryType geom_type() const override { return GeometryType::CSG; }
2,147,483,647✔
394

UNCOV
395
protected:
×
396
  //! Returns the beginning position of a parenthesis block (immediately before
397
  //! two surface tokens) in the RPN given a starting position at the end of
398
  //! that block (immediately after two surface tokens)
399
  //! \param start Starting position of the search
400
  //! \param rpn The rpn being searched
401
  static vector<int32_t>::iterator find_left_parenthesis(
402
    vector<int32_t>::iterator start, const vector<int32_t>& rpn);
403

404
private:
405
  Region region_;
406
};
407

408
//==============================================================================
409
//! Define an instance of a particular cell
410
//==============================================================================
411

412
//!  Stores information used to identify a unique cell in the model
413
struct CellInstance {
414
  //! Check for equality
415
  bool operator==(const CellInstance& other) const
3,900,232✔
416
  {
417
    return index_cell == other.index_cell && instance == other.instance;
3,900,232✔
418
  }
419

420
  int64_t index_cell;
421
  int64_t instance;
422
};
423

424
//! Structure necessary for inserting CellInstance into hashed STL data
425
//! structures
426
struct CellInstanceHash {
427
  std::size_t operator()(const CellInstance& k) const
3,900,882✔
428
  {
429
    return 4096 * k.index_cell + k.instance;
3,900,882✔
430
  }
431
};
432

433
//==============================================================================
434
// Non-member functions
435
//==============================================================================
436

437
void read_cells(pugi::xml_node node);
438

439
//!  Add cells to universes
440
void populate_universes();
441

442
} // namespace openmc
443
#endif // OPENMC_CELL_H
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