/
cell.h
405 lines (320 loc) · 14.4 KB
/
cell.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
#ifndef OPENMC_CELL_H
#define OPENMC_CELL_H
#include <cstdint>
#include <functional> // for hash
#include <limits>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include "hdf5.h"
#include "pugixml.hpp"
#include <gsl/gsl-lite.hpp>
#include "openmc/constants.h"
#include "openmc/memory.h" // for unique_ptr
#include "openmc/neighbor_list.h"
#include "openmc/position.h"
#include "openmc/surface.h"
#include "openmc/universe.h"
#include "openmc/vector.h"
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
enum class Fill { MATERIAL, UNIVERSE, LATTICE };
// TODO: Convert to enum
constexpr int32_t OP_LEFT_PAREN {std::numeric_limits<int32_t>::max()};
constexpr int32_t OP_RIGHT_PAREN {std::numeric_limits<int32_t>::max() - 1};
constexpr int32_t OP_COMPLEMENT {std::numeric_limits<int32_t>::max() - 2};
constexpr int32_t OP_INTERSECTION {std::numeric_limits<int32_t>::max() - 3};
constexpr int32_t OP_UNION {std::numeric_limits<int32_t>::max() - 4};
//==============================================================================
// Global variables
//==============================================================================
class Cell;
class GeometryState;
class ParentCell;
class CellInstance;
class Universe;
class UniversePartitioner;
namespace model {
extern std::unordered_map<int32_t, int32_t> cell_map;
extern vector<unique_ptr<Cell>> cells;
} // namespace model
//==============================================================================
class Region {
public:
//----------------------------------------------------------------------------
// Constructors
Region() {}
explicit Region(std::string region_spec, int32_t cell_id);
//----------------------------------------------------------------------------
// Methods
//! \brief Determine if a cell contains the particle at a given location.
//!
//! The bounds of the cell are determined by a logical expression involving
//! surface half-spaces. The expression used is given in infix notation
//!
//! The function is split into two cases, one for simple cells (those
//! involving only the intersection of half-spaces) and one for complex cells.
//! Both cases use short circuiting; however, in the case fo complex cells,
//! the complexity increases with the binary operators involved.
//! \param r The 3D Cartesian coordinate to check.
//! \param u A direction used to "break ties" the coordinates are very
//! close to a surface.
//! \param on_surface The signed index of a surface that the coordinate is
//! known to be on. This index takes precedence over surface sense
//! calculations.
bool contains(Position r, Direction u, int32_t on_surface) const;
//! Find the oncoming boundary of this cell.
std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface) const;
//! Get the BoundingBox for this cell.
BoundingBox bounding_box(int32_t cell_id) const;
//! Get the CSG expression as a string
std::string str() const;
//! Get a vector containing all the surfaces in the region expression
vector<int32_t> surfaces() const;
//----------------------------------------------------------------------------
// Accessors
//! Get Boolean of if the cell is simple or not
bool is_simple() const { return simple_; }
private:
//----------------------------------------------------------------------------
// Private Methods
//! Get a vector of the region expression in postfix notation
vector<int32_t> generate_postfix(int32_t cell_id) const;
//! Determine if a particle is inside the cell for a simple cell (only
//! intersection operators)
bool contains_simple(Position r, Direction u, int32_t on_surface) const;
//! Determine if a particle is inside the cell for a complex cell.
//!
//! Uses the comobination of half-spaces and binary operators to determine
//! if short circuiting can be used. Short cicuiting uses the relative and
//! absolute depth of parentheses in the expression.
bool contains_complex(Position r, Direction u, int32_t on_surface) const;
//! BoundingBox if the paritcle is in a simple cell.
BoundingBox bounding_box_simple() const;
//! BoundingBox if the particle is in a complex cell.
BoundingBox bounding_box_complex(vector<int32_t> postfix) const;
//! Enfource precedence: Parenthases, Complement, Intersection, Union
void add_precedence();
//! Add parenthesis to enforce precedence
std::vector<int32_t>::iterator add_parentheses(
std::vector<int32_t>::iterator start);
//! Remove complement operators from the expression
void remove_complement_ops();
//! Remove complement operators by using DeMorgan's laws
void apply_demorgan(
vector<int32_t>::iterator start, vector<int32_t>::iterator stop);
//----------------------------------------------------------------------------
// Private Data
//! Definition of spatial region as Boolean expression of half-spaces
// TODO: Should this be a vector of some other type
vector<int32_t> expression_;
bool simple_; //!< Does the region contain only intersections?
};
//==============================================================================
class Cell {
public:
//----------------------------------------------------------------------------
// Constructors, destructors, factory functions
explicit Cell(pugi::xml_node cell_node);
Cell() {};
virtual ~Cell() = default;
//----------------------------------------------------------------------------
// Methods
//! \brief Determine if a cell contains the particle at a given location.
//!
//! The bounds of the cell are detemined by a logical expression involving
//! surface half-spaces. At initialization, the expression was converted
//! to RPN notation.
//!
//! The function is split into two cases, one for simple cells (those
//! involving only the intersection of half-spaces) and one for complex cells.
//! Simple cells can be evaluated with short circuit evaluation, i.e., as soon
//! as we know that one half-space is not satisfied, we can exit. This
//! provides a performance benefit for the common case. In
//! contains_complex, we evaluate the RPN expression using a stack, similar to
//! how a RPN calculator would work.
//! \param r The 3D Cartesian coordinate to check.
//! \param u A direction used to "break ties" the coordinates are very
//! close to a surface.
//! \param on_surface The signed index of a surface that the coordinate is
//! known to be on. This index takes precedence over surface sense
//! calculations.
virtual bool contains(Position r, Direction u, int32_t on_surface) const = 0;
//! Find the oncoming boundary of this cell.
virtual std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, GeometryState* p) const = 0;
//! Write all information needed to reconstruct the cell to an HDF5 group.
//! \param group_id An HDF5 group id.
void to_hdf5(hid_t group_id) const;
virtual void to_hdf5_inner(hid_t group_id) const = 0;
//! Export physical properties to HDF5
//! \param[in] group HDF5 group to read from
void export_properties_hdf5(hid_t group) const;
//! Import physical properties from HDF5
//! \param[in] group HDF5 group to write to
void import_properties_hdf5(hid_t group);
//! Get the BoundingBox for this cell.
virtual BoundingBox bounding_box() const = 0;
//! Get a vector of surfaces in the cell
virtual vector<int32_t> surfaces() const { return vector<int32_t>(); }
//! Check if the cell region expression is simple
virtual bool is_simple() const { return true; }
//----------------------------------------------------------------------------
// Accessors
//! Get the temperature of a cell instance
//! \param[in] instance Instance index. If -1 is given, the temperature for
//! the first instance is returned.
//! \return Temperature in [K]
double temperature(int32_t instance = -1) const;
//! Set the temperature of a cell instance
//! \param[in] T Temperature in [K]
//! \param[in] instance Instance index. If -1 is given, the temperature for
//! all instances is set.
//! \param[in] set_contained If this cell is not filled with a material,
//! collect all contained cells with material fills and set their
//! temperatures.
void set_temperature(
double T, int32_t instance = -1, bool set_contained = false);
//! Set the rotation matrix of a cell instance
//! \param[in] rot The rotation matrix of length 3 or 9
void set_rotation(const vector<double>& rot);
//! Get the name of a cell
//! \return Cell name
const std::string& name() const { return name_; };
//! Set the temperature of a cell instance
//! \param[in] name Cell name
void set_name(const std::string& name) { name_ = name; };
//! Get all cell instances contained by this cell
//! \param[in] instance Instance of the cell for which to get contained cells
//! (default instance is zero)
//! \param[in] hint positional hint for determining the parent cells
//! \return Map with cell indexes as keys and
//! instances as values
std::unordered_map<int32_t, vector<int32_t>> get_contained_cells(
int32_t instance = 0, Position* hint = nullptr) const;
protected:
//! Determine the path to this cell instance in the geometry hierarchy
//! \param[in] instance of the cell to find parent cells for
//! \param[in] r position used to do a fast search for parent cells
//! \return parent cells
vector<ParentCell> find_parent_cells(
int32_t instance, const Position& r) const;
//! Determine the path to this cell instance in the geometry hierarchy
//! \param[in] instance of the cell to find parent cells for
//! \param[in] p particle used to do a fast search for parent cells
//! \return parent cells
vector<ParentCell> find_parent_cells(
int32_t instance, GeometryState& p) const;
//! Determine the path to this cell instance in the geometry hierarchy
//! \param[in] instance of the cell to find parent cells for
//! \return parent cells
vector<ParentCell> exhaustive_find_parent_cells(int32_t instance) const;
//! Inner function for retrieving contained cells
void get_contained_cells_inner(
std::unordered_map<int32_t, vector<int32_t>>& contained_cells,
vector<ParentCell>& parent_cells) const;
public:
//----------------------------------------------------------------------------
// Data members
int32_t id_; //!< Unique ID
std::string name_; //!< User-defined name
Fill type_; //!< Material, universe, or lattice
int32_t universe_; //!< Universe # this cell is in
int32_t fill_; //!< Universe # filling this cell
int32_t n_instances_ {0}; //!< Number of instances of this cell
GeometryType geom_type_; //!< Geometric representation type (CSG, DAGMC)
//! \brief Index corresponding to this cell in distribcell arrays
int distribcell_index_ {C_NONE};
//! \brief Material(s) within this cell.
//!
//! May be multiple materials for distribcell.
vector<int32_t> material_;
//! \brief Temperature(s) within this cell.
//!
//! The stored values are actually sqrt(k_Boltzmann * T) for each temperature
//! T. The units are sqrt(eV).
vector<double> sqrtkT_;
//! \brief Neighboring cells in the same universe.
NeighborList neighbors_;
Position translation_ {0, 0, 0}; //!< Translation vector for filled universe
//! \brief Rotational tranfsormation of the filled universe.
//
//! The vector is empty if there is no rotation. Otherwise, the first 9 values
//! give the rotation matrix in row-major order. When the user specifies
//! rotation angles about the x-, y- and z- axes in degrees, these values are
//! also present at the end of the vector, making it of length 12.
vector<double> rotation_;
vector<int32_t> offset_; //!< Distribcell offset table
};
struct CellInstanceItem {
int32_t index {-1}; //! Index into global cells array
int lattice_indx {-1}; //! Flat index value of the lattice cell
};
//==============================================================================
class CSGCell : public Cell {
public:
//----------------------------------------------------------------------------
// Constructors
CSGCell();
explicit CSGCell(pugi::xml_node cell_node);
//----------------------------------------------------------------------------
// Methods
vector<int32_t> surfaces() const override { return region_.surfaces(); }
std::pair<double, int32_t> distance(Position r, Direction u,
int32_t on_surface, GeometryState* p) const override
{
return region_.distance(r, u, on_surface);
}
bool contains(Position r, Direction u, int32_t on_surface) const override
{
return region_.contains(r, u, on_surface);
}
BoundingBox bounding_box() const override
{
return region_.bounding_box(id_);
}
void to_hdf5_inner(hid_t group_id) const override;
bool is_simple() const override { return region_.is_simple(); }
protected:
//! Returns the beginning position of a parenthesis block (immediately before
//! two surface tokens) in the RPN given a starting position at the end of
//! that block (immediately after two surface tokens)
//! \param start Starting position of the search
//! \param rpn The rpn being searched
static vector<int32_t>::iterator find_left_parenthesis(
vector<int32_t>::iterator start, const vector<int32_t>& rpn);
private:
Region region_;
};
//==============================================================================
//! Define an instance of a particular cell
//==============================================================================
//! Stores information used to identify a unique cell in the model
struct CellInstance {
//! Check for equality
bool operator==(const CellInstance& other) const
{
return index_cell == other.index_cell && instance == other.instance;
}
gsl::index index_cell;
gsl::index instance;
};
//! Structure necessary for inserting CellInstance into hashed STL data
//! structures
struct CellInstanceHash {
std::size_t operator()(const CellInstance& k) const
{
return 4096 * k.index_cell + k.instance;
}
};
//==============================================================================
// Non-member functions
//==============================================================================
void read_cells(pugi::xml_node node);
//! Add cells to universes
void populate_universes();
} // namespace openmc
#endif // OPENMC_CELL_H