From 3c1827c7dde907ccaae4acd1ecf90cab0f6fb289 Mon Sep 17 00:00:00 2001 From: shriwise Date: Thu, 17 May 2018 10:57:13 -0500 Subject: [PATCH] Abstracting cell class to support both CAD and CSG. --- include/openmc/cell.h | 48 +++++++++++++++++++++++++++++++++++++++++-- src/cad.cpp | 8 ++++---- src/cell.cpp | 31 ++++++++++++++++++++-------- 3 files changed, 72 insertions(+), 15 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index f819d41f44..9bbc004605 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -110,12 +110,49 @@ public: std::vector 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 + 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 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 diff --git a/src/cad.cpp b/src/cad.cpp index 65cab484f8..be008174e9 100644 --- a/src/cad.cpp +++ b/src/cad.cpp @@ -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); diff --git a/src/cell.cpp b/src/cell.cpp index fba9c8934b..6060edc7ad 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -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 -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::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 CADCell::distance(const double xyz[3], const double uvw[3], int32_t on_surface) const { + std::cout << "Cell Distance" << std::endl; + std::pair 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(); }