Abstracting cell class to support both CAD and CSG.

This commit is contained in:
shriwise 2018-05-17 10:57:13 -05:00 committed by pshriwise
parent df342211cc
commit 3c1827c7dd
3 changed files with 72 additions and 15 deletions

View file

@ -110,12 +110,49 @@ public:
std::vector<int32_t> offset_; //!< Distribcell offset table
Cell() {};
explicit Cell(pugi::xml_node cell_node);
explicit Cell();
//! 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 xyz[3] The 3D Cartesian coordinate to check.
//! @param uvw[3] 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(const double xyz[3], const double uvw[3], int32_t on_surface) const = 0;
virtual std::pair<double, int32_t>
distance(const double xyz[3], const double uvw[3], int32_t on_surface) const = 0;
//! Write all information needed to reconstruct the cell to an HDF5 group.
//! @param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const = 0;
};
class CSGCell : public Cell
{
public:
CSGCell();
explicit CSGCell(pugi::xml_node cell_node);
//! \brief Determine if a cell contains the particle at a given location.
//! Determine if a cell contains the particle at a given location.
//!
@ -147,6 +184,7 @@ public:
//! \param group_id An HDF5 group id.
void to_hdf5(hid_t group_id) const;
protected:
bool contains_simple(Position r, Direction u, int32_t on_surface) const;
bool contains_complex(Position r, Direction u, int32_t on_surface) const;
@ -158,6 +196,12 @@ class CADCell : public Cell
public:
moab::DagMC *dagmc_ptr;
explicit CADCell();
std::pair<double, int32_t> distance(const double xyz[3], const double uvw[3], int32_t on_surface) const;
bool contains(const double xyz[3], const double uvw[3], int32_t on_surface) const;
void to_hdf5(hid_t group_id) const;
};
#endif

View file

@ -23,12 +23,12 @@ void load_cad_geometry_c()
{
// set cell ids using global IDs
openmc::CADCell* c = new openmc::CADCell();
c->id = DAGMC->id_by_index(3, i);
c->id_ = DAGMC->id_by_index(3, i);
c->dagmc_ptr = DAGMC;
c->universe = cad_univ_id; // set to zero for now
c->material.push_back(40); // TEMPORARY
c->universe_ = cad_univ_id; // set to zero for now
c->material_.push_back(40); // TEMPORARY
openmc::cells_c.push_back(c);
openmc::cell_dict[c->id] = i;
openmc::cell_dict[c->id_] = i;
// Populate the Universe vector and dict
auto it = openmc::universe_dict.find(cad_univ_id);

View file

@ -213,9 +213,9 @@ Universe::to_hdf5(hid_t universes_group) const
// Cell implementation
//==============================================================================
Cell::Cell() {} // empty constructor
CSGCell::CSGCell() {} // empty constructor
Cell::Cell(pugi::xml_node cell_node)
CSGCell::CSGCell(pugi::xml_node cell_node)
{
if (check_for_node(cell_node, "id")) {
id_ = std::stoi(get_node_value(cell_node, "id"));
@ -399,7 +399,7 @@ Cell::Cell(pugi::xml_node cell_node)
//==============================================================================
bool
Cell::contains(Position r, Direction u, int32_t on_surface) const
CSGCell::contains(Position r, Direction u, int32_t on_surface) const
{
if (simple_) {
return contains_simple(r, u, on_surface);
@ -411,7 +411,7 @@ Cell::contains(Position r, Direction u, int32_t on_surface) const
//==============================================================================
std::pair<double, int32_t>
Cell::distance(Position r, Direction u, int32_t on_surface) const
CSGCell::distance(Position r, Direction u, int32_t on_surface) const
{
double min_dist {INFTY};
int32_t i_surf {std::numeric_limits<int32_t>::max()};
@ -440,7 +440,7 @@ Cell::distance(Position r, Direction u, int32_t on_surface) const
//==============================================================================
void
Cell::to_hdf5(hid_t cells_group) const
CSGCell::to_hdf5(hid_t cell_group) const
{
// Create a group for this cell.
std::stringstream group_name;
@ -519,7 +519,7 @@ Cell::to_hdf5(hid_t cells_group) const
//==============================================================================
bool
Cell::contains_simple(Position r, Direction u, int32_t on_surface) const
CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const
{
for (int32_t token : rpn_) {
if (token < OP_UNION) {
@ -543,7 +543,7 @@ Cell::contains_simple(Position r, Direction u, int32_t on_surface) const
//==============================================================================
bool
Cell::contains_complex(Position r, Direction u, int32_t on_surface) const
CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const
{
// Make a stack of booleans. We don't know how big it needs to be, but we do
// know that rpn.size() is an upper-bound.
@ -596,6 +596,19 @@ Cell::contains_complex(Position r, Direction u, int32_t on_surface) const
//==============================================================================
#ifdef CAD
CADCell::CADCell() : Cell{} {};
std::pair<double, int32_t> CADCell::distance(const double xyz[3], const double uvw[3], int32_t on_surface) const {
std::cout << "Cell Distance" << std::endl;
std::pair<double, int> result(1.0, 0);
return result;
}
bool CADCell::contains(const double xyz[3], const double uvw[3], int32_t on_surface) const {
std::cout << "Cell Contains" << std::endl;
return true;
}
void CADCell::to_hdf5(hid_t group_id) const { return; }
#endif
//==============================================================================
@ -612,7 +625,7 @@ read_cells(pugi::xml_node* node)
// Loop over XML cell elements and populate the array.
cells.reserve(n_cells);
for (pugi::xml_node cell_node: node->children("cell")) {
cells.push_back(new Cell(cell_node));
cells_c.push_back(new CSGCell(cell_node));
}
// Populate the Universe vector and map.
@ -764,7 +777,7 @@ extern "C" {
{
cells.reserve(cells.size() + n);
for (int32_t i = 0; i < n; i++) {
cells.push_back(new Cell());
cells_c.push_back(new CSGCell());
}
n_cells = cells.size();
}