From b68e7ca9021b05a563a1c4c37a3c38e98b931618 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 6 May 2019 16:16:05 -0500 Subject: [PATCH 01/40] Adding bounding box methods to all Surface types. --- include/openmc/surface.h | 27 +++++++++++++-------------- src/surface.cpp | 35 ++++++++++++++++++++++++++++++++++- 2 files changed, 47 insertions(+), 15 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 7b0ccefd1..78e02d994 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -31,6 +31,7 @@ extern "C" const int BC_PERIODIC; //============================================================================== class Surface; +struct BoundingBox; namespace model { extern std::vector> surfaces; @@ -105,6 +106,8 @@ public: //TODO: this probably needs to include i_periodic for PeriodicSurface virtual void to_hdf5(hid_t group_id) const = 0; + //! Get the BoundingBox for this surface. + virtual BoundingBox bounding_box() const = 0; }; class CSGSurface : public Surface @@ -126,9 +129,7 @@ protected: class DAGSurface : public Surface { public: - moab::DagMC* dagmc_ptr_; DAGSurface(); - int32_t dag_index_; double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; @@ -137,7 +138,8 @@ public: //! Get the bounding box of this surface. BoundingBox bounding_box() const; - void to_hdf5(hid_t group_id) const; + moab::DagMC* dagmc_ptr_; + int32_t dag_index_; }; #endif //============================================================================== @@ -165,9 +167,6 @@ public: //! boundary condition. virtual bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const = 0; - - //! Get the bounding box for this surface. - virtual BoundingBox bounding_box() const = 0; }; //============================================================================== @@ -269,7 +268,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - + BoundingBox bounding_box() const; double y0_, z0_, radius_; }; @@ -288,7 +287,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - + BoundingBox bounding_box() const; double x0_, z0_, radius_; }; @@ -307,7 +306,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - + BoundingBox bounding_box() const; double x0_, y0_, radius_; }; @@ -326,7 +325,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - + BoundingBox bounding_box() const; double x0_, y0_, z0_, radius_; }; @@ -345,7 +344,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - + BoundingBox bounding_box() const; double x0_, y0_, z0_, radius_sq_; }; @@ -364,7 +363,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - + BoundingBox bounding_box() const; double x0_, y0_, z0_, radius_sq_; }; @@ -383,7 +382,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - + BoundingBox bounding_box() const; double x0_, y0_, z0_, radius_sq_; }; @@ -401,7 +400,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - + BoundingBox bounding_box() const; // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_; }; diff --git a/src/surface.cpp b/src/surface.cpp index f364694a0..7b610031f 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -645,7 +645,6 @@ Direction SurfaceXCylinder::normal(Position r) const return axis_aligned_cylinder_normal<0, 1, 2>(r, y0_, z0_); } - void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-cylinder", false); @@ -653,6 +652,9 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceXCylinder::bounding_box() const { + return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_}; +} //============================================================================== // SurfaceYCylinder implementation //============================================================================== @@ -686,6 +688,10 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceYCylinder::bounding_box() const { + return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_, z0_ + radius_}; +} + //============================================================================== // SurfaceZCylinder implementation //============================================================================== @@ -719,6 +725,11 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceZCylinder::bounding_box() const { + return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY, INFTY}; +} + + //============================================================================== // SurfaceSphere implementation //============================================================================== @@ -787,6 +798,12 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceSphere::bounding_box() const { + return {x0_ - radius_, x0_ + radius_, + y0_ - radius_, y0_ + radius_, + z0_ - radius_, z0_ + radius_}; +} + //============================================================================== // Generic functions for x-, y-, and z-, cones //============================================================================== @@ -905,6 +922,10 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceXCone::bounding_box() const { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; +} + //============================================================================== // SurfaceYCone implementation //============================================================================== @@ -938,6 +959,10 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceYCone::bounding_box() const { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; +} + //============================================================================== // SurfaceZCone implementation //============================================================================== @@ -971,6 +996,10 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceZCone::bounding_box() const { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; +} + //============================================================================== // SurfaceQuadric implementation //============================================================================== @@ -1064,6 +1093,10 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceQuadric::bounding_box() const { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; +} + //============================================================================== void read_surfaces(pugi::xml_node node) From 446e78cb87a68cb451839bf6cd642dd1174138dd Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 20 Jun 2019 09:54:13 -0500 Subject: [PATCH 02/40] Updating bounding box definitions for all surfaces. --- include/openmc/cell.h | 5 ++ include/openmc/surface.h | 61 +++++++++++++++-------- src/cell.cpp | 7 +++ src/surface.cpp | 102 ++++++++++++++++++++++++++++++--------- 4 files changed, 131 insertions(+), 44 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index afe40fe43..8329bc510 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -114,6 +114,9 @@ public: //! \param group_id An HDF5 group id. virtual void to_hdf5(hid_t group_id) const = 0; + //! Get the BoundingBox for this cell. + virtual BoundingBox bounding_box() const = 0; + //---------------------------------------------------------------------------- // Accessors @@ -192,6 +195,8 @@ public: void to_hdf5(hid_t group_id) const; + BoundingBox bounding_box() 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; diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 78e02d994..3c5ece4f9 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -44,12 +44,33 @@ namespace model { struct BoundingBox { - double xmin; - double xmax; - double ymin; - double ymax; - double zmin; - double zmax; + double xmin = INFTY; + double xmax = -INFTY; + double ymin = INFTY; + double ymax = -INFTY; + double zmin = INFTY; + double zmax = -INFTY; + + // in-place update + inline void update(const BoundingBox& other) { + xmin = std::min(xmin, other.xmin); + xmax = std::max(xmax, other.xmax); + ymin = std::min(ymin, other.ymin); + ymax = std::max(ymax, other.ymax); + zmin = std::min(zmin, other.zmin); + zmax = std::max(zmax, other.zmax); + }; + + // in-place intersection + inline void intersect(const BoundingBox& other) { + xmin = std::max(xmin, other.xmin); + xmax = std::min(xmax, other.xmax); + ymin = std::max(ymin, other.ymin); + ymax = std::min(ymax, other.ymax); + zmin = std::max(zmin, other.zmin); + zmax = std::min(zmax, other.zmax); + } + }; //============================================================================== @@ -107,7 +128,7 @@ public: virtual void to_hdf5(hid_t group_id) const = 0; //! Get the BoundingBox for this surface. - virtual BoundingBox bounding_box() const = 0; + virtual BoundingBox bounding_box(bool pos_side) const = 0; }; class CSGSurface : public Surface @@ -136,7 +157,7 @@ public: Direction normal(Position r) const; Direction reflect(Position r, Direction u) const; //! Get the bounding box of this surface. - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; moab::DagMC* dagmc_ptr_; int32_t dag_index_; @@ -185,7 +206,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double x0_; }; @@ -206,7 +227,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double y0_; }; @@ -227,7 +248,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double z0_; }; @@ -248,7 +269,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double A_, B_, C_, D_; }; @@ -268,7 +289,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double y0_, z0_, radius_; }; @@ -287,7 +308,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double x0_, z0_, radius_; }; @@ -306,7 +327,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, radius_; }; @@ -325,7 +346,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, z0_, radius_; }; @@ -344,7 +365,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, z0_, radius_sq_; }; @@ -363,7 +384,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, z0_, radius_sq_; }; @@ -382,7 +403,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, z0_, radius_sq_; }; @@ -400,7 +421,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_; }; diff --git a/src/cell.cpp b/src/cell.cpp index 04273eded..9afae2d60 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -567,6 +567,13 @@ CSGCell::to_hdf5(hid_t cell_group) const close_group(group); } +BoundingBox CSGCell::bounding_box() const { + BoundingBox bbox; + for (int32_t token : rpn_) { + bbox.update(model::surfaces[abs(token)-1]->bounding_box(token > 0)); + } +} + //============================================================================== bool diff --git a/src/surface.cpp b/src/surface.cpp index 7b610031f..e9d2f964e 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -271,7 +271,7 @@ Direction DAGSurface::reflect(Position r, Direction u) const return simulation::last_dir; } -BoundingBox DAGSurface::bounding_box() const +BoundingBox DAGSurface::bounding_box(bool pos_side) const { moab::ErrorCode rval; moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_); @@ -366,9 +366,13 @@ bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other, } BoundingBox -SurfaceXPlane::bounding_box() const +SurfaceXPlane::bounding_box(bool pos_side) const { - return {x0_, x0_, -INFTY, INFTY, -INFTY, INFTY}; + if (pos_side) { + return {x0_, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + } else { + return {-INFTY, x0_, -INFTY, INFTY, -INFTY, INFTY}; + } } //============================================================================== @@ -428,9 +432,13 @@ bool SurfaceYPlane::periodic_translate(const PeriodicSurface* other, } BoundingBox -SurfaceYPlane::bounding_box() const +SurfaceYPlane::bounding_box(bool pos_side) const { - return {-INFTY, INFTY, y0_, y0_, -INFTY, INFTY}; + if (pos_side) { + return {-INFTY, INFTY, y0_, INFTY, -INFTY, INFTY}; + } else { + return {-INFTY, INFTY, -INFTY, y0_, -INFTY, INFTY}; + } } //============================================================================== @@ -474,9 +482,13 @@ bool SurfaceZPlane::periodic_translate(const PeriodicSurface* other, } BoundingBox -SurfaceZPlane::bounding_box() const +SurfaceZPlane::bounding_box(bool pos_side) const { - return {-INFTY, INFTY, -INFTY, INFTY, z0_, z0_}; + if (pos_side) { + return {-INFTY, INFTY, -INFTY, INFTY, z0_, INFTY}; + } else { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, z0_}; + } } //============================================================================== @@ -539,9 +551,34 @@ bool SurfacePlane::periodic_translate(const PeriodicSurface* other, Position& r, } BoundingBox -SurfacePlane::bounding_box() const +SurfacePlane::bounding_box(bool pos_side) const { - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + BoundingBox bbox = {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + if (A_ == 0.0 && B_ == 0.0) { + double val = D_ / C_; + if (pos_side) { + bbox.zmin = val; + } else { + bbox.zmax = val; + } + } else if (A_ == 0.0 && C_ == 0.0) { + double val = D_ / B_; + if (pos_side) { + bbox.ymin = val; + } else { + bbox.ymax = val; + } + } else if (B_ == 0.0 && C_ == 0.0) { + double val = D_ / A_; + if (pos_side) { + bbox.xmin = val; + } else { + bbox.xmax = val; + } + } + + return bbox; + } //============================================================================== @@ -652,8 +689,12 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceXCylinder::bounding_box() const { - return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_}; +BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const { + if (pos_side) { + return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_}; + } else { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + } } //============================================================================== // SurfaceYCylinder implementation @@ -688,8 +729,12 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceYCylinder::bounding_box() const { - return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_, z0_ + radius_}; +BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const { + if (pos_side) { + return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_, z0_ + radius_}; + } else { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + } } //============================================================================== @@ -725,8 +770,12 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceZCylinder::bounding_box() const { - return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY, INFTY}; +BoundingBox SurfaceZCylinder::bounding_box(bool pos_side) const { + if (pos_side) { + return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY, INFTY}; + } else { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + } } @@ -798,10 +847,14 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceSphere::bounding_box() const { - return {x0_ - radius_, x0_ + radius_, - y0_ - radius_, y0_ + radius_, - z0_ - radius_, z0_ + radius_}; +BoundingBox SurfaceSphere::bounding_box(bool pos_side) const { + if (pos_side) { + return {x0_ - radius_, x0_ + radius_, + y0_ - radius_, y0_ + radius_, + z0_ - radius_, z0_ + radius_}; + } else { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + } } //============================================================================== @@ -922,7 +975,7 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceXCone::bounding_box() const { +BoundingBox SurfaceXCone::bounding_box(bool pos_side) const { return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; } @@ -959,7 +1012,7 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceYCone::bounding_box() const { +BoundingBox SurfaceYCone::bounding_box(bool pos_side) const { return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; } @@ -996,7 +1049,7 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceZCone::bounding_box() const { +BoundingBox SurfaceZCone::bounding_box(bool pos_side) const { return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; } @@ -1093,7 +1146,7 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceQuadric::bounding_box() const { +BoundingBox SurfaceQuadric::bounding_box(bool pos_side) const { return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; } @@ -1194,7 +1247,8 @@ void read_surfaces(pugi::xml_node node) } // See if this surface makes part of the global bounding box. - BoundingBox bb = surf->bounding_box(); + BoundingBox bb = surf->bounding_box(true); + bb.intersect(surf->bounding_box(false)); if (bb.xmin > -INFTY && bb.xmin < xmin) { xmin = bb.xmin; i_xmin = i_surf; From 470e26c78337585bb2b5e703ba4a464d3177865e Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 21 Jun 2019 00:27:18 -0500 Subject: [PATCH 03/40] Adding bbox defs for cells. --- include/openmc/cell.h | 2 ++ include/openmc/surface.h | 7 +++++++ src/cell.cpp | 10 ++++++++++ src/surface.cpp | 1 + 4 files changed, 20 insertions(+) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 8329bc510..b4848bb50 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -217,6 +217,8 @@ public: std::pair distance(Position r, Direction u, int32_t on_surface) const; + BoundingBox bounding_box() const; + void to_hdf5(hid_t group_id) const; }; #endif diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 3c5ece4f9..f218f7ac7 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -139,6 +139,8 @@ public: void to_hdf5(hid_t group_id) const; + virtual BoundingBox bounding_box(bool pos_side) const = 0; + protected: virtual void to_hdf5_inner(hid_t group_id) const = 0; }; @@ -159,6 +161,8 @@ public: //! Get the bounding box of this surface. BoundingBox bounding_box(bool pos_side) const; + void to_hdf5(hid_t group_id) const; + moab::DagMC* dagmc_ptr_; int32_t dag_index_; }; @@ -188,6 +192,9 @@ public: //! boundary condition. virtual bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const = 0; + + //! Get the bounding box for this surface. + virtual BoundingBox bounding_box(bool pos_side) const = 0; }; //============================================================================== diff --git a/src/cell.cpp b/src/cell.cpp index 9afae2d60..ebf466bd6 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -698,6 +698,16 @@ bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const void DAGCell::to_hdf5(hid_t group_id) const { return; } +BoundingBox DAGCell::bounding_box() const +{ + moab::ErrorCode rval; + moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_); + double min[3], max[3]; + rval = dagmc_ptr_->getobb(vol, min, max); + MB_CHK_ERR_CONT(rval); + return {min[0], max[0], min[1], max[1], min[2], max[2]}; +} + #endif //============================================================================== diff --git a/src/surface.cpp b/src/surface.cpp index e9d2f964e..ae8a5b6c4 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -282,6 +282,7 @@ BoundingBox DAGSurface::bounding_box(bool pos_side) const } void DAGSurface::to_hdf5(hid_t group_id) const {} + #endif //============================================================================== // PeriodicSurface implementation From 08239a98d276e897f0bdffec6826727531817841 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 21 Jun 2019 14:51:06 -0500 Subject: [PATCH 04/40] Adding a bounding box method to the universe class. --- include/openmc/cell.h | 2 ++ src/cell.cpp | 22 ++++++++++++++++++++-- 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index b4848bb50..a6432a304 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -65,6 +65,8 @@ public: //! \param group_id An HDF5 group id. void to_hdf5(hid_t group_id) const; + BoundingBox bounding_box() const; + std::unique_ptr partitioner_; }; diff --git a/src/cell.cpp b/src/cell.cpp index ebf466bd6..107724c61 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -208,6 +208,19 @@ Universe::to_hdf5(hid_t universes_group) const close_group(group); } +BoundingBox Universe::bounding_box() const { + BoundingBox bbox; + if (cells_.size() == 0) { + bbox = {-INFTY, INFTY, -INFTY, -INFTY, INFTY}; + } else { + for (const auto& cell : cells_) { + auto& c = model::cells[cell]; + bbox.update(c->bounding_box()); + } + } + return bbox; +} + //============================================================================== // Cell implementation //============================================================================== @@ -569,9 +582,14 @@ CSGCell::to_hdf5(hid_t cell_group) const BoundingBox CSGCell::bounding_box() const { BoundingBox bbox; - for (int32_t token : rpn_) { - bbox.update(model::surfaces[abs(token)-1]->bounding_box(token > 0)); + if (rpn_.size() == 0) { + bbox = {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + } else { + for (int32_t token : rpn_) { + bbox.update(model::surfaces[abs(token)-1]->bounding_box(token > 0)); + } } + return bbox; } //============================================================================== From ecbf7596c0f84f889e98969351070f4b2717803d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 25 Jun 2019 12:33:22 -0500 Subject: [PATCH 05/40] Corrections to surface bounding boxes based on sign. Corrections to cell/universe bounding boxes and capi exposure. --- include/openmc/capi.h | 1 + include/openmc/surface.h | 16 +++++++++------ openmc/capi/core.py | 17 ++++++++++++++++ src/cell.cpp | 4 ++-- src/geometry.cpp | 43 ++++++++++++++++++++++++++++++++++++++++ src/surface.cpp | 8 ++++---- 6 files changed, 77 insertions(+), 12 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 148a28c88..ad2adae24 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -29,6 +29,7 @@ extern "C" { int openmc_filter_set_id(int32_t index, int32_t id); int openmc_finalize(); int openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance); + int openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double* urc); int openmc_fission_bank(void** ptr, int64_t* n); int openmc_get_cell_index(int32_t id, int32_t* index); int openmc_get_filter_index(int32_t id, int32_t* index); diff --git a/include/openmc/surface.h b/include/openmc/surface.h index f218f7ac7..69c11476c 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -44,12 +44,12 @@ namespace model { struct BoundingBox { - double xmin = INFTY; - double xmax = -INFTY; - double ymin = INFTY; - double ymax = -INFTY; - double zmin = INFTY; - double zmax = -INFTY; + double xmin = -INFTY; + double xmax = INFTY; + double ymin = -INFTY; + double ymax = INFTY; + double zmin = -INFTY; + double zmax = INFTY; // in-place update inline void update(const BoundingBox& other) { @@ -71,8 +71,12 @@ struct BoundingBox zmax = std::min(zmax, other.zmax); } + }; + + + //============================================================================== //! A geometry primitive used to define regions of 3D space. //============================================================================== diff --git a/openmc/capi/core.py b/openmc/capi/core.py index aae17e06a..5744e123c 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -74,6 +74,23 @@ _dll.openmc_statepoint_write.argtypes = [c_char_p, POINTER(c_bool)] _dll.openmc_statepoint_write.restype = c_int _dll.openmc_statepoint_write.errcheck = _error_handler +_dll.openmc_bounding_box.argtypes = [c_char_p, c_int, POINTER(c_double), + POINTER(c_double)] +_dll.openmc_bounding_box.restype = c_int +_dll.openmc_bounding_box.errcheck = _error_handler + +def bounding_box(geom_type, geom_id): + + geomt = c_char_p(geom_type.encode()) + llc = np.zeros((3,), dtype=float) + urc = np.zeros((3,), dtype=float) + _dll.openmc_bounding_box(geomt, + geom_id, + llc.ctypes.data_as(POINTER(c_double)), + urc.ctypes.data_as(POINTER(c_double))) + + return llc, urc + def calculate_volumes(): """Run stochastic volume calculation""" diff --git a/src/cell.cpp b/src/cell.cpp index 107724c61..d646660fb 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -209,7 +209,7 @@ Universe::to_hdf5(hid_t universes_group) const } BoundingBox Universe::bounding_box() const { - BoundingBox bbox; + BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY}; if (cells_.size() == 0) { bbox = {-INFTY, INFTY, -INFTY, -INFTY, INFTY}; } else { @@ -586,7 +586,7 @@ BoundingBox CSGCell::bounding_box() const { bbox = {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; } else { for (int32_t token : rpn_) { - bbox.update(model::surfaces[abs(token)-1]->bounding_box(token > 0)); + bbox.intersect(model::surfaces[abs(token)-1]->bounding_box(token > 0)); } } return bbox; diff --git a/src/geometry.cpp b/src/geometry.cpp index 0f89a0ecd..ed75444ec 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -9,6 +9,7 @@ #include "openmc/lattice.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/string_utils.h" #include "openmc/surface.h" @@ -475,4 +476,46 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance) return 0; } +extern "C" int +openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double* urc) { + + BoundingBox bbox; + + std::string gtype(geom_type); + to_lower(gtype); + + std::cout << gtype << std::endl; + + if (gtype == "universe") { + // negative ids only apply to surfaces + if (id < 0) { return OPENMC_E_GEOMETRY; } + const auto& u = model::universes[model::universe_map[id]]; + bbox = u->bounding_box(); + } else if (gtype == "cell") { + // negative ids only apply to surfaces + if (id < 0) { return OPENMC_E_GEOMETRY; } + const auto& c = model::cells[model::cell_map[id]]; + bbox = c->bounding_box(); + } else if (gtype == "surface") { + const auto& s = model::surfaces[model::surface_map[abs(id)]]; + bbox = s->bounding_box(id > 0); + } else { + std::stringstream msg; + msg << "Geometry type: " << gtype << " is invalid."; + return OPENMC_E_GEOMETRY; + } + + // set lower left corner values + llc[0] = bbox.xmin; + llc[1] = bbox.ymin; + llc[2] = bbox.zmin; + + // set upper right corner values + urc[0] = bbox.xmax; + urc[1] = bbox.ymax; + urc[2] = bbox.zmax; + + return 0; +} + } // namespace openmc diff --git a/src/surface.cpp b/src/surface.cpp index ae8a5b6c4..9eaccef3d 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -691,7 +691,7 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const } BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const { - if (pos_side) { + if (!pos_side) { return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_}; } else { return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; @@ -731,7 +731,7 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const } BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const { - if (pos_side) { + if (!pos_side) { return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_, z0_ + radius_}; } else { return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; @@ -772,7 +772,7 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const } BoundingBox SurfaceZCylinder::bounding_box(bool pos_side) const { - if (pos_side) { + if (!pos_side) { return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY, INFTY}; } else { return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; @@ -849,7 +849,7 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const } BoundingBox SurfaceSphere::bounding_box(bool pos_side) const { - if (pos_side) { + if (!pos_side) { return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_}; From 2ac505814e22019644459eedac70aa9eb24e7d08 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 9 Jul 2019 11:11:25 -0500 Subject: [PATCH 06/40] Adding test for universe bounding box. --- tests/unit_tests/test_capi.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index af86a054d..aa938f8b2 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -469,3 +469,18 @@ def test_position(capi_init): pos[2] = 3.3 assert tuple(pos) == (1.3, 2.3, 3.3) + +def test_bounding_box(capi_init): + + expected_llc = (-0.63, -0.63, -np.inf) + expected_urc = (0.63, 0.63, np.inf) + + llc, urc = openmc.capi.bounding_box("Universe", 0) + + assert llc[0] == expected_llc[0] + assert llc[1] == expected_llc[1] + assert llc[2] < 1.0E10 + + assert urc[0] == expected_urc[0] + assert urc[1] == expected_urc[1] + assert urc[2] > 1.E10 From a10570245520dee4d1aafcfcc42a356ffa7c8f91 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 9 Jul 2019 16:23:24 -0500 Subject: [PATCH 07/40] Updates to CAPI bounding box methods. Addition of global bounding box method. Adding tests for PWR pincell model. --- include/openmc/capi.h | 1 + openmc/capi/core.py | 14 +++++++- src/geometry.cpp | 24 ++++++++++--- tests/unit_tests/test_capi.py | 67 ++++++++++++++++++++++++++++++----- 4 files changed, 92 insertions(+), 14 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index ad2adae24..11b8a6ef6 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -30,6 +30,7 @@ extern "C" { int openmc_finalize(); int openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance); int openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double* urc); + int openmc_global_bounding_box(double* llc, double* urc); int openmc_fission_bank(void** ptr, int64_t* n); int openmc_get_cell_index(int32_t id, int32_t* index); int openmc_get_filter_index(int32_t id, int32_t* index); diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 5744e123c..ed78b38b6 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -73,11 +73,23 @@ _dll.openmc_simulation_finalize.errcheck = _error_handler _dll.openmc_statepoint_write.argtypes = [c_char_p, POINTER(c_bool)] _dll.openmc_statepoint_write.restype = c_int _dll.openmc_statepoint_write.errcheck = _error_handler - _dll.openmc_bounding_box.argtypes = [c_char_p, c_int, POINTER(c_double), POINTER(c_double)] _dll.openmc_bounding_box.restype = c_int _dll.openmc_bounding_box.errcheck = _error_handler +_dll.openmc_global_bounding_box.argtypes = [POINTER(c_double), + POINTER(c_double)] +_dll.openmc_global_bounding_box.restype = c_int +_dll.openmc_global_bounding_box.errcheck = _error_handler + +def global_bounding_box(): + + llc = np.zeros((3,), dtype=float) + urc = np.zeros((3,), dtype=float) + _dll.openmc_global_bounding_box(llc.ctypes.data_as(POINTER(c_double)), + urc.ctypes.data_as(POINTER(c_double))) + + return llc, urc def bounding_box(geom_type, geom_id): diff --git a/src/geometry.cpp b/src/geometry.cpp index ed75444ec..43620644f 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -484,19 +484,18 @@ openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double std::string gtype(geom_type); to_lower(gtype); - std::cout << gtype << std::endl; - if (gtype == "universe") { // negative ids only apply to surfaces - if (id < 0) { return OPENMC_E_GEOMETRY; } + if (id <= 0) { return OPENMC_E_GEOMETRY; } const auto& u = model::universes[model::universe_map[id]]; bbox = u->bounding_box(); } else if (gtype == "cell") { // negative ids only apply to surfaces - if (id < 0) { return OPENMC_E_GEOMETRY; } + if (id <= 0) { return OPENMC_E_GEOMETRY; } const auto& c = model::cells[model::cell_map[id]]; bbox = c->bounding_box(); } else if (gtype == "surface") { + if (id == 0) { return OPENMC_E_GEOMETRY; } const auto& s = model::surfaces[model::surface_map[abs(id)]]; bbox = s->bounding_box(id > 0); } else { @@ -518,4 +517,21 @@ openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double return 0; } +extern "C" int openmc_global_bounding_box(double* llc, double* urc) { + auto bbox = model::universes[model::root_universe]->bounding_box(); + + // set lower left corner values + llc[0] = bbox.xmin; + llc[1] = bbox.ymin; + llc[2] = bbox.zmin; + + // set upper right corner values + urc[0] = bbox.xmax; + urc[1] = bbox.ymax; + urc[2] = bbox.zmax; + + return 0; +} + + } // namespace openmc diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index aa938f8b2..d611c20c7 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -1,5 +1,6 @@ from collections.abc import Mapping import os +import sys import numpy as np import pytest @@ -472,15 +473,63 @@ def test_position(capi_init): def test_bounding_box(capi_init): - expected_llc = (-0.63, -0.63, -np.inf) - expected_urc = (0.63, 0.63, np.inf) + inf = sys.float_info.max - llc, urc = openmc.capi.bounding_box("Universe", 0) + expected_llc = (-inf, -0.63, -inf) + expected_urc = (inf, inf, inf) - assert llc[0] == expected_llc[0] - assert llc[1] == expected_llc[1] - assert llc[2] < 1.0E10 + llc, urc = openmc.capi.bounding_box("Surface", 5) - assert urc[0] == expected_urc[0] - assert urc[1] == expected_urc[1] - assert urc[2] > 1.E10 + print(llc) + print(urc) + + assert tuple(llc) == expected_llc + assert tuple(urc) == expected_urc + + + expected_llc = (-inf, -inf, -inf) + expected_urc = (inf, -0.63, inf) + + llc, urc = openmc.capi.bounding_box("Surface", -5) + + assert tuple(llc) == expected_llc + assert tuple(urc) == expected_urc + + expected_llc = (-0.39218, -0.39218, -inf) + expected_urc = (0.39218, 0.39218, inf) + + llc, urc = openmc.capi.bounding_box("Cell", 1) + + assert tuple(llc) == expected_llc + assert tuple(urc) == expected_urc + + expected_llc = (-0.45720, -0.45720, -inf) + expected_urc = (0.45720, 0.45720, inf) + + llc, urc = openmc.capi.bounding_box("Cell", 2) + + assert tuple(llc) == expected_llc + assert tuple(urc) == expected_urc + + # make sure that proper assertions are raised + with pytest.raises(openmc.exceptions.GeometryError): + openmc.capi.bounding_box("Cell", -1) + + with pytest.raises(openmc.exceptions.GeometryError): + openmc.capi.bounding_box("Surface", 0) + + with pytest.raises(openmc.exceptions.GeometryError): + openmc.capi.bounding_box("Region", 1) + + +def test_global_bounding_box(capi_init): + + inf = sys.float_info.max + + expected_llc = (-0.63, -0.63, -inf) + expected_urc = (0.63, 0.63, inf) + + llc, urc = openmc.capi.global_bounding_box() + + assert tuple(llc) == expected_llc + assert tuple(urc) == expected_urc From 1d075cefdd2f295b917c274a9ee165cf6f7d3123 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 10 Jul 2019 07:53:44 -0500 Subject: [PATCH 08/40] Adding name property to CAPI cell/material classes. --- include/openmc/capi.h | 4 ++++ openmc/capi/cell.py | 19 ++++++++++++++++++- openmc/capi/material.py | 19 ++++++++++++++++++- src/cell.cpp | 26 ++++++++++++++++++++++++++ src/material.cpp | 25 +++++++++++++++++++++++++ 5 files changed, 91 insertions(+), 2 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 11b8a6ef6..d81e1e722 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -14,6 +14,8 @@ extern "C" { int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n); int openmc_cell_get_id(int32_t index, int32_t* id); int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T); + int openmc_cell_get_name(int32_t index, const char*& name); + int openmc_cell_set_name(int32_t index, const char* name); int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices); int openmc_cell_set_id(int32_t index, int32_t id); int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance); @@ -58,6 +60,8 @@ extern "C" { int openmc_material_set_density(int32_t index, double density, const char* units); int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density); int openmc_material_set_id(int32_t index, int32_t id); + int openmc_material_get_name(int32_t index, const char*& name); + int openmc_material_set_name(int32_t index, const char* name); int openmc_material_set_volume(int32_t index, double volume); int openmc_material_filter_get_bins(int32_t index, const int32_t** bins, size_t* n); int openmc_material_filter_set_bins(int32_t index, size_t n, const int32_t* bins); diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 959ab08fc..5f38e823e 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -1,5 +1,5 @@ from collections.abc import Mapping, Iterable -from ctypes import c_int, c_int32, c_double, c_char_p, POINTER +from ctypes import byref, c_int, c_int32, c_double, c_char_p, POINTER from weakref import WeakValueDictionary import numpy as np @@ -28,6 +28,12 @@ _dll.openmc_cell_get_temperature.argtypes = [ c_int32, POINTER(c_int32), POINTER(c_double)] _dll.openmc_cell_get_temperature.restype = c_int _dll.openmc_cell_get_temperature.errcheck = _error_handler +_dll.openmc_cell_get_name.argtypes = [c_int32, POINTER(c_char_p)] +_dll.openmc_cell_get_name.restype = c_int +_dll.openmc_cell_get_name.errcheck = _error_handler +_dll.openmc_cell_set_name.argtypes = [c_int32, c_char_p] +_dll.openmc_cell_set_name.restype = c_int +_dll.openmc_cell_set_name.errcheck = _error_handler _dll.openmc_cell_set_fill.argtypes = [ c_int32, c_int, c_int32, POINTER(c_int32)] _dll.openmc_cell_set_fill.restype = c_int @@ -102,6 +108,17 @@ class Cell(_FortranObjectWithID): def id(self, cell_id): _dll.openmc_cell_set_id(self._index, cell_id) + @property + def name(self): + name = c_char_p() + _dll.openmc_cell_get_name(self._index, byref(name)) + return name.value.decode() + + @name.setter + def name(self, name): + name_ptr = c_char_p(name.encode()) + _dll.openmc_cell_set_name(self._index, name_ptr) + @property def fill(self): fill_type = c_int() diff --git a/openmc/capi/material.py b/openmc/capi/material.py index 43959e244..3f98e9f67 100644 --- a/openmc/capi/material.py +++ b/openmc/capi/material.py @@ -1,5 +1,5 @@ from collections.abc import Mapping -from ctypes import c_int, c_int32, c_double, c_char_p, POINTER, c_size_t +from ctypes import byref, c_int, c_int32, c_double, c_char_p, POINTER, c_size_t from weakref import WeakValueDictionary import numpy as np @@ -48,6 +48,12 @@ _dll.openmc_material_set_densities.errcheck = _error_handler _dll.openmc_material_set_id.argtypes = [c_int32, c_int32] _dll.openmc_material_set_id.restype = c_int _dll.openmc_material_set_id.errcheck = _error_handler +_dll.openmc_cell_get_name.argtypes = [c_int32, POINTER(c_char_p)] +_dll.openmc_cell_get_name.restype = c_int +_dll.openmc_cell_get_name.errcheck = _error_handler +_dll.openmc_cell_set_name.argtypes = [c_int32, c_char_p] +_dll.openmc_cell_set_name.restype = c_int +_dll.openmc_cell_set_name.errcheck = _error_handler _dll.openmc_material_set_volume.argtypes = [c_int32, c_double] _dll.openmc_material_set_volume.restype = c_int _dll.openmc_material_set_volume.errcheck = _error_handler @@ -124,6 +130,17 @@ class Material(_FortranObjectWithID): def id(self, mat_id): _dll.openmc_material_set_id(self._index, mat_id) + @property + def name(self): + name = c_char_p() + _dll.openmc_material_get_name(self._index, byref(name)) + return name.value.decode() + + @name.setter + def name(self, name): + name_ptr = c_char_p(name.encode()) + _dll.openmc_material_set_name(self._index, name_ptr) + @property def volume(self): volume = c_double() diff --git a/src/cell.cpp b/src/cell.cpp index d646660fb..5c1c90db5 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -1018,6 +1018,32 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) return 0; } +extern "C" int +openmc_cell_get_name(int32_t index, const char*& name) { + if (index < 0 || index >= model::cells.size()) { + strcpy(openmc_err_msg, "Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + name = model::cells[index]->name_.c_str(); + + return 0; +} + +extern "C" int +openmc_cell_set_name(int32_t index, const char* name) { + if (index < 0 || index >= model::cells.size()) { + strcpy(openmc_err_msg, "Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + std::string name_str(name); + model::cells[index]->name_ = name_str; + + return 0; +} + + //! Return the index in the cells array of a cell with a given ID extern "C" int openmc_get_cell_index(int32_t id, int32_t* index) diff --git a/src/material.cpp b/src/material.cpp index 0a6bebf17..fb29ba2b7 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1381,6 +1381,31 @@ openmc_material_set_id(int32_t index, int32_t id) return 0; } +extern "C" int +openmc_material_get_name(int32_t index, const char*& name) { + if (index < 0 || index >= model::materials.size()) { + strcpy(openmc_err_msg, "Index in materials array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + name = model::materials[index]->name_.c_str(); + + return 0; +} + +extern "C" int +openmc_material_set_name(int32_t index, const char* name) { + if (index < 0 || index >= model::materials.size()) { + strcpy(openmc_err_msg, "Index in materials array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + std::string name_str(name); + model::materials[index]->name_ = name_str; + + return 0; +} + extern "C" int openmc_material_set_volume(int32_t index, double volume) { From 542ce3f15e9a85e23103dd69ee71457e173a11f8 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 Jul 2019 21:52:09 -0500 Subject: [PATCH 09/40] Adding material/cell name tests. --- tests/unit_tests/test_capi.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index d611c20c7..a5e0e6a20 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -74,7 +74,7 @@ def test_cell(capi_init): assert isinstance(cell.fill, openmc.capi.Material) cell.fill = openmc.capi.materials[1] assert str(cell) == 'Cell[0]' - + assert cell.name == "Fuel" def test_cell_temperature(capi_init): cell = openmc.capi.cells[1] @@ -124,6 +124,7 @@ def test_material(capi_init): m.set_density(0.1, 'g/cm3') assert m.density == pytest.approx(0.1) + assert m.name == "Hot borated water" def test_material_add_nuclide(capi_init): m = openmc.capi.materials[3] From e9ac00736748fbfa09cef3d396ef749047cc8166 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 12 Jul 2019 15:52:37 -0500 Subject: [PATCH 10/40] Removing complement from internal RPN. --- include/openmc/cell.h | 2 ++ include/openmc/surface.h | 3 -- src/cell.cpp | 71 ++++++++++++++++++++++++++++++++-------- 3 files changed, 59 insertions(+), 17 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index a6432a304..37db57bcf 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -202,6 +202,8 @@ public: 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; + BoundingBox bounding_box_simple() const; + BoundingBox bounding_box_complex() const; }; //============================================================================== diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 69c11476c..2b6c5cd0c 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -74,9 +74,6 @@ struct BoundingBox }; - - - //============================================================================== //! A geometry primitive used to define regions of 3D space. //============================================================================== diff --git a/src/cell.cpp b/src/cell.cpp index 5c1c90db5..aa937dae7 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -114,13 +114,19 @@ generate_rpn(int32_t cell_id, std::vector infix) std::vector rpn; std::vector stack; + bool in_complement = false; + for (int32_t token : infix) { if (token < OP_UNION) { // If token is not an operator, add it to output + if (in_complement) { token *= -1; } rpn.push_back(token); - } else if (token < OP_RIGHT_PAREN) { // Regular operators union, intersection, complement + if (in_complement) { + if (token == OP_UNION) token = OP_INTERSECTION; + else if (token == OP_INTERSECTION) token = OP_UNION; + } while (stack.size() > 0) { int32_t op = stack.back(); @@ -133,14 +139,19 @@ generate_rpn(int32_t cell_id, std::vector infix) // is less than that of op, move op to the output queue and push the // token on to the stack. Note that only complement is // right-associative. - rpn.push_back(op); + if (op == OP_COMPLEMENT) { in_complement = false; } + else { rpn.push_back(op); } stack.pop_back(); } else { break; } } - stack.push_back(token); + if (token == OP_COMPLEMENT) { + in_complement = true; + } else { + stack.push_back(token); + } } else if (token == OP_LEFT_PAREN) { // If the token is a left parenthesis, push it onto the stack @@ -158,8 +169,9 @@ generate_rpn(int32_t cell_id, std::vector infix) << cell_id; fatal_error(err_msg); } - - rpn.push_back(stack.back()); + int32_t op = stack.back(); + if (op == OP_COMPLEMENT) { in_complement = false; } + else { rpn.push_back(op); } stack.pop_back(); } @@ -183,6 +195,10 @@ generate_rpn(int32_t cell_id, std::vector infix) stack.pop_back(); } + for (int32_t val : rpn) { + std::cout << val << std::endl; + } + return rpn; } @@ -580,18 +596,47 @@ CSGCell::to_hdf5(hid_t cell_group) const close_group(group); } -BoundingBox CSGCell::bounding_box() const { +BoundingBox CSGCell::bounding_box_simple() const { BoundingBox bbox; - if (rpn_.size() == 0) { - bbox = {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; - } else { - for (int32_t token : rpn_) { - bbox.intersect(model::surfaces[abs(token)-1]->bounding_box(token > 0)); - } + for (int32_t token : rpn_) { + bbox.intersect(model::surfaces[abs(token)-1]->bounding_box(token > 0)); } return bbox; } +BoundingBox CSGCell::bounding_box_complex() const { + + std::vector stack(rpn_.size()); + int i_stack = -1; + + for (int32_t token : rpn_) { + if (token == OP_UNION) { + stack[i_stack - 1].update(stack[i_stack]); + i_stack--; + } else if (token == OP_INTERSECTION) { + stack[i_stack - 1].intersect(stack[i_stack]); + i_stack--; + } else { + i_stack++; + stack[i_stack] = model::surfaces[abs(token)-1]->bounding_box(token > 0); + } + } + return stack[i_stack]; +} + +BoundingBox CSGCell::bounding_box() const { + if (rpn_.size() == 0) { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + } else { + if (simple_) { + return bounding_box_simple(); + } + else { + return bounding_box_complex(); + } + } +} + //============================================================================== bool @@ -634,8 +679,6 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const } else if (token == OP_INTERSECTION) { stack[i_stack-1] = stack[i_stack-1] && stack[i_stack]; i_stack --; - } else if (token == OP_COMPLEMENT) { - stack[i_stack] = !stack[i_stack]; } else { // If the token is not an operator, evaluate the sense of particle with // respect to the surface and see if the token matches the sense. If the From d223046b30760ed6847b9363c1bc51655c7147a1 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 15 Jul 2019 12:50:21 -0500 Subject: [PATCH 11/40] Adding complement operator back into internal cell RPN. Updating complex_cell bounding box algorithm to parse existing RPN. --- include/openmc/cell.h | 2 +- src/cell.cpp | 93 ++++++++++++++++++++++++++----------------- 2 files changed, 58 insertions(+), 37 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 37db57bcf..aadd038d3 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -203,7 +203,7 @@ 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; BoundingBox bounding_box_simple() const; - BoundingBox bounding_box_complex() const; + BoundingBox bounding_box_complex(std::vector rpn) const; }; //============================================================================== diff --git a/src/cell.cpp b/src/cell.cpp index aa937dae7..586ccf220 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -114,19 +114,12 @@ generate_rpn(int32_t cell_id, std::vector infix) std::vector rpn; std::vector stack; - bool in_complement = false; - for (int32_t token : infix) { if (token < OP_UNION) { // If token is not an operator, add it to output - if (in_complement) { token *= -1; } rpn.push_back(token); } else if (token < OP_RIGHT_PAREN) { // Regular operators union, intersection, complement - if (in_complement) { - if (token == OP_UNION) token = OP_INTERSECTION; - else if (token == OP_INTERSECTION) token = OP_UNION; - } while (stack.size() > 0) { int32_t op = stack.back(); @@ -139,19 +132,14 @@ generate_rpn(int32_t cell_id, std::vector infix) // is less than that of op, move op to the output queue and push the // token on to the stack. Note that only complement is // right-associative. - if (op == OP_COMPLEMENT) { in_complement = false; } - else { rpn.push_back(op); } + rpn.push_back(op); stack.pop_back(); } else { break; } } - if (token == OP_COMPLEMENT) { - in_complement = true; - } else { - stack.push_back(token); - } + stack.push_back(token); } else if (token == OP_LEFT_PAREN) { // If the token is a left parenthesis, push it onto the stack @@ -169,9 +157,7 @@ generate_rpn(int32_t cell_id, std::vector infix) << cell_id; fatal_error(err_msg); } - int32_t op = stack.back(); - if (op == OP_COMPLEMENT) { in_complement = false; } - else { rpn.push_back(op); } + rpn.push_back(stack.back()); stack.pop_back(); } @@ -195,10 +181,6 @@ generate_rpn(int32_t cell_id, std::vector infix) stack.pop_back(); } - for (int32_t val : rpn) { - std::cout << val << std::endl; - } - return rpn; } @@ -604,24 +586,61 @@ BoundingBox CSGCell::bounding_box_simple() const { return bbox; } -BoundingBox CSGCell::bounding_box_complex() const { +BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { - std::vector stack(rpn_.size()); - int i_stack = -1; + std::reverse(rpn.begin(), rpn.end()); - for (int32_t token : rpn_) { - if (token == OP_UNION) { - stack[i_stack - 1].update(stack[i_stack]); - i_stack--; - } else if (token == OP_INTERSECTION) { - stack[i_stack - 1].intersect(stack[i_stack]); - i_stack--; - } else { - i_stack++; - stack[i_stack] = model::surfaces[abs(token)-1]->bounding_box(token > 0); + BoundingBox current = model::surfaces[abs(rpn.back()) - 1]->bounding_box(rpn.back() > 0); + rpn.pop_back(); + + while (rpn.size()) { + int32_t one = rpn.back(); rpn.pop_back(); + int32_t two = rpn.back(); rpn.pop_back(); + + assert(one < OP_UNION); + + if (two >= OP_UNION) { + if (two == OP_UNION) { + current.update(model::surfaces[abs(one)-1]->bounding_box(one > 0)); + } else if (two == OP_INTERSECTION) { + current.intersect(model::surfaces[abs(one)-1]->bounding_box(one > 0)); + } + } else { // two surfaces in a row, create sub-rpn for parenthesis + std::vector subrpn; + subrpn.push_back(one); + subrpn.push_back(two); + int32_t sone = one; + int32_t stwo = two; + // add until last two tokens are operators + while ((subrpn.back() < OP_UNION) || (*(subrpn.rbegin() + 1) < OP_UNION)) { + subrpn.push_back(rpn.back()); + rpn.pop_back(); + } + + // handle complement case + if (subrpn.back() == OP_COMPLEMENT) { + subrpn.pop_back(); + for (auto& token : subrpn) { + if (token < OP_UNION) { token *= -1; } + else if (token == OP_UNION) { token = OP_INTERSECTION; } + else if (token == OP_INTERSECTION) { token = OP_UNION; } + } + subrpn.push_back(rpn.back()); + rpn.pop_back(); + } + // get bbox for the subrpn + int32_t op = subrpn.back(); + subrpn.pop_back(); + BoundingBox sub_box = bounding_box_complex(subrpn); + if (op == OP_UNION) { + current.update(sub_box); + } else if (op == OP_INTERSECTION) { + current.intersect(sub_box); + } } } - return stack[i_stack]; + + return current; } BoundingBox CSGCell::bounding_box() const { @@ -632,7 +651,7 @@ BoundingBox CSGCell::bounding_box() const { return bounding_box_simple(); } else { - return bounding_box_complex(); + return bounding_box_complex(rpn_); } } } @@ -679,6 +698,8 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const } else if (token == OP_INTERSECTION) { stack[i_stack-1] = stack[i_stack-1] && stack[i_stack]; i_stack --; + } else if (token == OP_COMPLEMENT) { + stack[i_stack] = !stack[i_stack]; } else { // If the token is not an operator, evaluate the sense of particle with // respect to the surface and see if the token matches the sense. If the From 77d40effb15aad4ee9d39883ca7b9bac8db75c0b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 15 Jul 2019 14:35:02 -0500 Subject: [PATCH 12/40] Adding a test for complex cell bounding boxes. --- openmc/capi/cell.py | 8 +++++++ tests/regression_tests/complex_cell/test.py | 25 +++++++++++++++++++++ 2 files changed, 33 insertions(+) diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 5f38e823e..1388d9732 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -183,6 +183,14 @@ class Cell(_FortranObjectWithID): _dll.openmc_cell_set_temperature(self._index, T, instance) + @property + def bounding_box(self): + llc = np.zeros((3,), dtype=float) + urc = np.zeros((3,), dtype=float) + _dll.openmc_bounding_box(b'Cell', self.id, + llc.ctypes.data_as(POINTER(c_double)), + urc.ctypes.data_as(POINTER(c_double))) + return llc, urc class _CellMapping(Mapping): def __getitem__(self, key): diff --git a/tests/regression_tests/complex_cell/test.py b/tests/regression_tests/complex_cell/test.py index 77cbd6cb7..816902b4c 100755 --- a/tests/regression_tests/complex_cell/test.py +++ b/tests/regression_tests/complex_cell/test.py @@ -1,6 +1,31 @@ from tests.testing_harness import TestHarness +import sys + +import openmc.capi def test_complex_cell(): harness = TestHarness('statepoint.10.h5') harness.main() + +def test_complex_cell_capi(): + # initialize + openmc.capi.init([]) + + inf = sys.float_info.max + + expected_boxes = { 1 : (( -4., -4., -inf), ( 4., 4., inf)), + 2 : (( -7., -7., -inf), ( 7., 7., inf)), + 3 : ((-10., -10., -inf), (10., 10., inf)), + 4 : ((-10., -10., -inf), (10., 10., inf)) } + + for cell_id, cell in openmc.capi.cells.items(): + cell_box = cell.bounding_box + + assert tuple(cell_box[0]) == expected_boxes[cell_id][0] + assert tuple(cell_box[1]) == expected_boxes[cell_id][1] + + cell_box = openmc.capi.bounding_box("Cell", cell_id) + + assert tuple(cell_box[0]) == expected_boxes[cell_id][0] + assert tuple(cell_box[1]) == expected_boxes[cell_id][1] From 9fab83c70c257ab02c1059412d3652ca97bb6b36 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 15 Jul 2019 14:35:20 -0500 Subject: [PATCH 13/40] Style fix in capi core. --- openmc/capi/core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index ed78b38b6..a3ae6a219 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -87,7 +87,7 @@ def global_bounding_box(): llc = np.zeros((3,), dtype=float) urc = np.zeros((3,), dtype=float) _dll.openmc_global_bounding_box(llc.ctypes.data_as(POINTER(c_double)), - urc.ctypes.data_as(POINTER(c_double))) + urc.ctypes.data_as(POINTER(c_double))) return llc, urc From 1d02c874660dea68c0ef248881f684ea8b239451 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 15 Jul 2019 15:02:07 -0500 Subject: [PATCH 14/40] Some self-review. --- include/openmc/surface.h | 13 +++++++++++-- openmc/capi/core.py | 3 ++- src/cell.cpp | 21 +++++++++++---------- 3 files changed, 24 insertions(+), 13 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 2b6c5cd0c..6a688c240 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -31,6 +31,7 @@ extern "C" const int BC_PERIODIC; //============================================================================== class Surface; + struct BoundingBox; namespace model { @@ -51,7 +52,7 @@ struct BoundingBox double zmin = -INFTY; double zmax = INFTY; - // in-place update + // in-place update to include another bounding box inline void update(const BoundingBox& other) { xmin = std::min(xmin, other.xmin); xmax = std::max(xmax, other.xmax); @@ -61,7 +62,7 @@ struct BoundingBox zmax = std::max(zmax, other.zmax); }; - // in-place intersection + // in-place intersection with another bounding box inline void intersect(const BoundingBox& other) { xmin = std::max(xmin, other.xmin); xmax = std::min(xmax, other.xmax); @@ -298,6 +299,7 @@ public: Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + double y0_, z0_, radius_; }; @@ -317,6 +319,7 @@ public: Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + double x0_, z0_, radius_; }; @@ -336,6 +339,7 @@ public: Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + double x0_, y0_, radius_; }; @@ -355,6 +359,7 @@ public: Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + double x0_, y0_, z0_, radius_; }; @@ -374,6 +379,7 @@ public: Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + double x0_, y0_, z0_, radius_sq_; }; @@ -393,6 +399,7 @@ public: Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + double x0_, y0_, z0_, radius_sq_; }; @@ -412,6 +419,7 @@ public: Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + double x0_, y0_, z0_, radius_sq_; }; @@ -430,6 +438,7 @@ public: Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; BoundingBox bounding_box(bool pos_side) const; + // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_; }; diff --git a/openmc/capi/core.py b/openmc/capi/core.py index a3ae6a219..d37e0de75 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -73,7 +73,8 @@ _dll.openmc_simulation_finalize.errcheck = _error_handler _dll.openmc_statepoint_write.argtypes = [c_char_p, POINTER(c_bool)] _dll.openmc_statepoint_write.restype = c_int _dll.openmc_statepoint_write.errcheck = _error_handler -_dll.openmc_bounding_box.argtypes = [c_char_p, c_int, POINTER(c_double), +_dll.openmc_bounding_box.argtypes = [c_char_p, c_int, + POINTER(c_double), POINTER(c_double)] _dll.openmc_bounding_box.restype = c_int _dll.openmc_bounding_box.errcheck = _error_handler diff --git a/src/cell.cpp b/src/cell.cpp index 586ccf220..a9c40ac24 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -594,9 +594,11 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { rpn.pop_back(); while (rpn.size()) { + // move through the rpn in twos int32_t one = rpn.back(); rpn.pop_back(); int32_t two = rpn.back(); rpn.pop_back(); + // the first token should always be a surface assert(one < OP_UNION); if (two >= OP_UNION) { @@ -605,19 +607,20 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { } else if (two == OP_INTERSECTION) { current.intersect(model::surfaces[abs(one)-1]->bounding_box(one > 0)); } - } else { // two surfaces in a row, create sub-rpn for parenthesis + } else { + // two surfaces in a row, create sub-rpn for region in parenthesis std::vector subrpn; subrpn.push_back(one); subrpn.push_back(two); int32_t sone = one; int32_t stwo = two; - // add until last two tokens are operators + // add until last two tokens in the sub-rpn are operators while ((subrpn.back() < OP_UNION) || (*(subrpn.rbegin() + 1) < OP_UNION)) { subrpn.push_back(rpn.back()); rpn.pop_back(); } - // handle complement case + // handle complement case using De Morgan's laws if (subrpn.back() == OP_COMPLEMENT) { subrpn.pop_back(); for (auto& token : subrpn) { @@ -628,10 +631,12 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { subrpn.push_back(rpn.back()); rpn.pop_back(); } - // get bbox for the subrpn - int32_t op = subrpn.back(); - subrpn.pop_back(); + // save the last operator, tells us how to combine this region + // with our current bounding box + int32_t op = subrpn.back(); subrpn.pop_back(); + // get bounding box for the subrpn BoundingBox sub_box = bounding_box_complex(subrpn); + // combine the sub-rpn bounding box with our current cell box if (op == OP_UNION) { current.update(sub_box); } else if (op == OP_INTERSECTION) { @@ -644,16 +649,12 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { } BoundingBox CSGCell::bounding_box() const { - if (rpn_.size() == 0) { - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; - } else { if (simple_) { return bounding_box_simple(); } else { return bounding_box_complex(rpn_); } - } } //============================================================================== From 3e0dc41ea14e8402ec6f3a16330f57002dd91e0f Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 15 Jul 2019 20:05:47 -0500 Subject: [PATCH 15/40] Adding documentation to CAPI bounding box function. --- openmc/capi/core.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index d37e0de75..686dd7862 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -84,7 +84,7 @@ _dll.openmc_global_bounding_box.restype = c_int _dll.openmc_global_bounding_box.errcheck = _error_handler def global_bounding_box(): - + """Calculate a global bounding box for the model""" llc = np.zeros((3,), dtype=float) urc = np.zeros((3,), dtype=float) _dll.openmc_global_bounding_box(llc.ctypes.data_as(POINTER(c_double)), @@ -93,7 +93,15 @@ def global_bounding_box(): return llc, urc def bounding_box(geom_type, geom_id): + """Get a bounding box for a geometric object + Parameters + ---------- + geom_type : str + Type of geometry object. One of ('surface', 'cell', 'universe') + geom_id : int + Id of the object. Can be positive or negative for surfaces. + """ geomt = c_char_p(geom_type.encode()) llc = np.zeros((3,), dtype=float) urc = np.zeros((3,), dtype=float) From b373fd5931dca10f223feee5b5e7b3454c12f85e Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 16 Jul 2019 10:11:58 -0500 Subject: [PATCH 16/40] Finalizing the capi after the complex cell test. --- tests/regression_tests/complex_cell/test.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/regression_tests/complex_cell/test.py b/tests/regression_tests/complex_cell/test.py index 816902b4c..8074c0e05 100755 --- a/tests/regression_tests/complex_cell/test.py +++ b/tests/regression_tests/complex_cell/test.py @@ -29,3 +29,5 @@ def test_complex_cell_capi(): assert tuple(cell_box[0]) == expected_boxes[cell_id][0] assert tuple(cell_box[1]) == expected_boxes[cell_id][1] + + openmc.capi.finalize() From 7409ad4632b219e20d26a9dad9207c19a86dfccf Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 16 Jul 2019 10:12:24 -0500 Subject: [PATCH 17/40] Some doc strings for new capi functions. --- src/cell.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/cell.cpp b/src/cell.cpp index a9c40ac24..1242c650f 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -1083,6 +1083,7 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) return 0; } +//! Get the name of a cell extern "C" int openmc_cell_get_name(int32_t index, const char*& name) { if (index < 0 || index >= model::cells.size()) { @@ -1095,6 +1096,7 @@ openmc_cell_get_name(int32_t index, const char*& name) { return 0; } +//! Set the name of a cell extern "C" int openmc_cell_set_name(int32_t index, const char* name) { if (index < 0 || index >= model::cells.size()) { From d79a5edb8a4bedfa2aadeed995b5c7d5339a05d8 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 16 Jul 2019 16:47:29 -0500 Subject: [PATCH 18/40] Small formatting fixes. --- tests/unit_tests/test_capi.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index a5e0e6a20..1c1da98a3 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -123,7 +123,6 @@ def test_material(capi_init): m.set_density(0.1, 'g/cm3') assert m.density == pytest.approx(0.1) - assert m.name == "Hot borated water" def test_material_add_nuclide(capi_init): @@ -473,7 +472,6 @@ def test_position(capi_init): assert tuple(pos) == (1.3, 2.3, 3.3) def test_bounding_box(capi_init): - inf = sys.float_info.max expected_llc = (-inf, -0.63, -inf) @@ -524,7 +522,6 @@ def test_bounding_box(capi_init): def test_global_bounding_box(capi_init): - inf = sys.float_info.max expected_llc = (-0.63, -0.63, -inf) From 08e78a3215e99f8ce6a0b8f36dc7c286f1a8badf Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 16 Jul 2019 17:09:44 -0500 Subject: [PATCH 19/40] Removing unecessary forward declaration. --- include/openmc/surface.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 6a688c240..e78c49079 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -32,8 +32,6 @@ extern "C" const int BC_PERIODIC; class Surface; -struct BoundingBox; - namespace model { extern std::vector> surfaces; extern std::unordered_map surface_map; From df379471dd2b6bf5f83ceb79c023e9aa595a7325 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 12:19:13 -0500 Subject: [PATCH 20/40] Adding Cell accessors for name. --- include/openmc/cell.h | 8 ++++++++ src/cell.cpp | 19 ++++++++++++++++--- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index aadd038d3..0b219baa1 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -134,6 +134,14 @@ public: //! all instances is set. void set_temperature(double T, int32_t instance = -1); + //! Get the name of a cell + //! \return Cell name + std::string name() const; + + //! Set the temperature of a cell instance + //! \param[in] name Cell name + void set_name(const std::string& name); + //---------------------------------------------------------------------------- // Data members diff --git a/src/cell.cpp b/src/cell.cpp index 1242c650f..c339a54c1 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -252,6 +252,18 @@ Cell::set_temperature(double T, int32_t instance) } } +std::string +Cell::name() const +{ + return name_; +} + +void +Cell::set_name(const std::string& name) +{ + name_ = name; +} + //============================================================================== // CSGCell implementation //============================================================================== @@ -1085,13 +1097,14 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) //! Get the name of a cell extern "C" int -openmc_cell_get_name(int32_t index, const char*& name) { +openmc_cell_get_name(int32_t index, char*& name) { if (index < 0 || index >= model::cells.size()) { strcpy(openmc_err_msg, "Index in cells array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - name = model::cells[index]->name_.c_str(); + auto name_str = model::cells[index]->name_; + strcpy(name, name_str.c_str()); return 0; } @@ -1105,7 +1118,7 @@ openmc_cell_set_name(int32_t index, const char* name) { } std::string name_str(name); - model::cells[index]->name_ = name_str; + model::cells[index]->set_name(name_str); return 0; } From 54df6518887d2eabfe8954fa1951826aed87d598 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 12:25:53 -0500 Subject: [PATCH 21/40] Adding Material accessors for name. --- include/openmc/cell.h | 4 ++-- include/openmc/material.h | 7 +++++++ src/cell.cpp | 12 ------------ src/material.cpp | 7 ++++--- 4 files changed, 13 insertions(+), 17 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 0b219baa1..3a0c81646 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -136,11 +136,11 @@ public: //! Get the name of a cell //! \return Cell name - std::string 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); + void set_name(const std::string& name) { name_ = name; }; //---------------------------------------------------------------------------- // Data members diff --git a/include/openmc/material.h b/include/openmc/material.h index 02c5380f3..aa0a4bad0 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -92,6 +92,13 @@ public: //! \return Density in [g/cm^3] double density_gpcc() const { return density_gpcc_; } + //! Get name + //! \return Material name + std::string name() const { return name_; } + + //! Set name + void set_name(const std::string& name) { name_ = name; } + //! Set total density of the material // //! \param[in] density Density value diff --git a/src/cell.cpp b/src/cell.cpp index c339a54c1..9d5b2af9d 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -252,18 +252,6 @@ Cell::set_temperature(double T, int32_t instance) } } -std::string -Cell::name() const -{ - return name_; -} - -void -Cell::set_name(const std::string& name) -{ - name_ = name; -} - //============================================================================== // CSGCell implementation //============================================================================== diff --git a/src/material.cpp b/src/material.cpp index fb29ba2b7..f14bfc8f2 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1382,13 +1382,14 @@ openmc_material_set_id(int32_t index, int32_t id) } extern "C" int -openmc_material_get_name(int32_t index, const char*& name) { +openmc_material_get_name(int32_t index, char*& name) { if (index < 0 || index >= model::materials.size()) { strcpy(openmc_err_msg, "Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - name = model::materials[index]->name_.c_str(); + auto name_str = model::materials[index]->name(); + strcpy(name, name_str.c_str()); return 0; } @@ -1401,7 +1402,7 @@ openmc_material_set_name(int32_t index, const char* name) { } std::string name_str(name); - model::materials[index]->name_ = name_str; + model::materials[index]->set_name(name); return 0; } From 520b07a131ce41962ba965a3d90085231fb626f6 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 12:29:30 -0500 Subject: [PATCH 22/40] Using C-types only in capi function definitions. --- src/cell.cpp | 4 ++-- src/material.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 9d5b2af9d..7753e8333 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -1085,14 +1085,14 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) //! Get the name of a cell extern "C" int -openmc_cell_get_name(int32_t index, char*& name) { +openmc_cell_get_name(int32_t index, char** name) { if (index < 0 || index >= model::cells.size()) { strcpy(openmc_err_msg, "Index in cells array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } auto name_str = model::cells[index]->name_; - strcpy(name, name_str.c_str()); + strcpy(*name, name_str.c_str()); return 0; } diff --git a/src/material.cpp b/src/material.cpp index f14bfc8f2..356c2bf21 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1382,14 +1382,14 @@ openmc_material_set_id(int32_t index, int32_t id) } extern "C" int -openmc_material_get_name(int32_t index, char*& name) { +openmc_material_get_name(int32_t index, char** name) { if (index < 0 || index >= model::materials.size()) { strcpy(openmc_err_msg, "Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } auto name_str = model::materials[index]->name(); - strcpy(name, name_str.c_str()); + strcpy(*name, name_str.c_str()); return 0; } From 6d971884139a5ca40a0874a68471a89ce170691c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 12:57:23 -0500 Subject: [PATCH 23/40] Creating a default virtual definition for a surfce bounding box and removing redundant definitions of surfaces with unbounded boxes. --- include/openmc/surface.h | 13 +-------- src/surface.cpp | 57 ---------------------------------------- 2 files changed, 1 insertion(+), 69 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index e78c49079..10b4a4c27 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -128,7 +128,7 @@ public: virtual void to_hdf5(hid_t group_id) const = 0; //! Get the BoundingBox for this surface. - virtual BoundingBox bounding_box(bool pos_side) const = 0; + virtual BoundingBox bounding_box(bool pos_side) const { return {}; } }; class CSGSurface : public Surface @@ -139,8 +139,6 @@ public: void to_hdf5(hid_t group_id) const; - virtual BoundingBox bounding_box(bool pos_side) const = 0; - protected: virtual void to_hdf5_inner(hid_t group_id) const = 0; }; @@ -158,8 +156,6 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; Direction reflect(Position r, Direction u) const; - //! Get the bounding box of this surface. - BoundingBox bounding_box(bool pos_side) const; void to_hdf5(hid_t group_id) const; @@ -193,8 +189,6 @@ public: virtual bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const = 0; - //! Get the bounding box for this surface. - virtual BoundingBox bounding_box(bool pos_side) const = 0; }; //============================================================================== @@ -276,7 +270,6 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; - BoundingBox bounding_box(bool pos_side) const; double A_, B_, C_, D_; }; @@ -376,7 +369,6 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, z0_, radius_sq_; }; @@ -396,7 +388,6 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, z0_, radius_sq_; }; @@ -416,7 +407,6 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, z0_, radius_sq_; }; @@ -435,7 +425,6 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - BoundingBox bounding_box(bool pos_side) const; // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_; diff --git a/src/surface.cpp b/src/surface.cpp index 9eaccef3d..b53590e08 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -271,16 +271,6 @@ Direction DAGSurface::reflect(Position r, Direction u) const return simulation::last_dir; } -BoundingBox DAGSurface::bounding_box(bool pos_side) const -{ - moab::ErrorCode rval; - moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_); - double min[3], max[3]; - rval = dagmc_ptr_->getobb(surf, min, max); - MB_CHK_ERR_CONT(rval); - return {min[0], max[0], min[1], max[1], min[2], max[2]}; -} - void DAGSurface::to_hdf5(hid_t group_id) const {} #endif @@ -551,37 +541,6 @@ bool SurfacePlane::periodic_translate(const PeriodicSurface* other, Position& r, return false; } -BoundingBox -SurfacePlane::bounding_box(bool pos_side) const -{ - BoundingBox bbox = {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; - if (A_ == 0.0 && B_ == 0.0) { - double val = D_ / C_; - if (pos_side) { - bbox.zmin = val; - } else { - bbox.zmax = val; - } - } else if (A_ == 0.0 && C_ == 0.0) { - double val = D_ / B_; - if (pos_side) { - bbox.ymin = val; - } else { - bbox.ymax = val; - } - } else if (B_ == 0.0 && C_ == 0.0) { - double val = D_ / A_; - if (pos_side) { - bbox.xmin = val; - } else { - bbox.xmax = val; - } - } - - return bbox; - -} - //============================================================================== // Generic functions for x-, y-, and z-, cylinders //============================================================================== @@ -976,10 +935,6 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceXCone::bounding_box(bool pos_side) const { - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; -} - //============================================================================== // SurfaceYCone implementation //============================================================================== @@ -1013,10 +968,6 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceYCone::bounding_box(bool pos_side) const { - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; -} - //============================================================================== // SurfaceZCone implementation //============================================================================== @@ -1050,10 +1001,6 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceZCone::bounding_box(bool pos_side) const { - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; -} - //============================================================================== // SurfaceQuadric implementation //============================================================================== @@ -1147,10 +1094,6 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -BoundingBox SurfaceQuadric::bounding_box(bool pos_side) const { - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; -} - //============================================================================== void read_surfaces(pugi::xml_node node) From 07c9083446b95fcebcd44d849cfdb808f724ed47 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 13:02:08 -0500 Subject: [PATCH 24/40] Using default constructor when returning unbounded boxes. --- src/cell.cpp | 4 ++-- src/surface.cpp | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 7753e8333..1dff5caf9 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -207,9 +207,9 @@ Universe::to_hdf5(hid_t universes_group) const } BoundingBox Universe::bounding_box() const { - BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY}; + BoundingBox bbox; if (cells_.size() == 0) { - bbox = {-INFTY, INFTY, -INFTY, -INFTY, INFTY}; + return {}; } else { for (const auto& cell : cells_) { auto& c = model::cells[cell]; diff --git a/src/surface.cpp b/src/surface.cpp index b53590e08..45cd025ba 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -653,7 +653,7 @@ BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const { if (!pos_side) { return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_}; } else { - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + return {}; } } //============================================================================== @@ -693,7 +693,7 @@ BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const { if (!pos_side) { return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_, z0_ + radius_}; } else { - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + return {}; } } @@ -734,7 +734,7 @@ BoundingBox SurfaceZCylinder::bounding_box(bool pos_side) const { if (!pos_side) { return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY, INFTY}; } else { - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + return {}; } } @@ -813,7 +813,7 @@ BoundingBox SurfaceSphere::bounding_box(bool pos_side) const { y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_}; } else { - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + return {}; } } From 81bf0a67daf332d51b1479a60bcc8077388300f2 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 13:15:26 -0500 Subject: [PATCH 25/40] Implementing intersection and union operators for BoundingBox. --- include/openmc/surface.h | 34 +++++++++++++++++++++++----------- src/cell.cpp | 12 ++++++------ src/surface.cpp | 3 +-- 3 files changed, 30 insertions(+), 19 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 10b4a4c27..a7d107ee3 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -50,26 +50,38 @@ struct BoundingBox double zmin = -INFTY; double zmax = INFTY; - // in-place update to include another bounding box - inline void update(const BoundingBox& other) { - xmin = std::min(xmin, other.xmin); - xmax = std::max(xmax, other.xmax); - ymin = std::min(ymin, other.ymin); - ymax = std::max(ymax, other.ymax); - zmin = std::min(zmin, other.zmin); - zmax = std::max(zmax, other.zmax); - }; - // in-place intersection with another bounding box - inline void intersect(const BoundingBox& other) { + inline BoundingBox operator &(const BoundingBox& other) { + BoundingBox result = *this; + return result &= other; + } + + inline BoundingBox operator |(const BoundingBox& other) { + BoundingBox result = *this; + return result |= other; + } + + // intersect operator + inline BoundingBox& operator &=(const BoundingBox& other) { xmin = std::max(xmin, other.xmin); xmax = std::min(xmax, other.xmax); ymin = std::max(ymin, other.ymin); ymax = std::min(ymax, other.ymax); zmin = std::max(zmin, other.zmin); zmax = std::min(zmax, other.zmax); + return *this; } + // union operator + inline BoundingBox& operator |=(const BoundingBox& other) { + xmin = std::min(xmin, other.xmin); + xmax = std::max(xmax, other.xmax); + ymin = std::min(ymin, other.ymin); + ymax = std::max(ymax, other.ymax); + zmin = std::min(zmin, other.zmin); + zmax = std::max(zmax, other.zmax); + return *this; + } }; diff --git a/src/cell.cpp b/src/cell.cpp index 1dff5caf9..fc9a25e2e 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -213,7 +213,7 @@ BoundingBox Universe::bounding_box() const { } else { for (const auto& cell : cells_) { auto& c = model::cells[cell]; - bbox.update(c->bounding_box()); + bbox |= c->bounding_box(); } } return bbox; @@ -581,7 +581,7 @@ CSGCell::to_hdf5(hid_t cell_group) const BoundingBox CSGCell::bounding_box_simple() const { BoundingBox bbox; for (int32_t token : rpn_) { - bbox.intersect(model::surfaces[abs(token)-1]->bounding_box(token > 0)); + bbox &= model::surfaces[abs(token)-1]->bounding_box(token > 0); } return bbox; } @@ -603,9 +603,9 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { if (two >= OP_UNION) { if (two == OP_UNION) { - current.update(model::surfaces[abs(one)-1]->bounding_box(one > 0)); + current |= model::surfaces[abs(one)-1]->bounding_box(one > 0); } else if (two == OP_INTERSECTION) { - current.intersect(model::surfaces[abs(one)-1]->bounding_box(one > 0)); + current &= model::surfaces[abs(one)-1]->bounding_box(one > 0); } } else { // two surfaces in a row, create sub-rpn for region in parenthesis @@ -638,9 +638,9 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { BoundingBox sub_box = bounding_box_complex(subrpn); // combine the sub-rpn bounding box with our current cell box if (op == OP_UNION) { - current.update(sub_box); + current |= sub_box; } else if (op == OP_INTERSECTION) { - current.intersect(sub_box); + current &= sub_box; } } } diff --git a/src/surface.cpp b/src/surface.cpp index 45cd025ba..14d077cff 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1191,8 +1191,7 @@ void read_surfaces(pugi::xml_node node) } // See if this surface makes part of the global bounding box. - BoundingBox bb = surf->bounding_box(true); - bb.intersect(surf->bounding_box(false)); + auto bb = surf->bounding_box(true) & surf->bounding_box(false); if (bb.xmin > -INFTY && bb.xmin < xmin) { xmin = bb.xmin; i_xmin = i_surf; From c1d8f06ac1ced6a9b40c6971538c4148dea6311b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 14:49:45 -0500 Subject: [PATCH 26/40] Restoring initial universe bounding box. --- src/cell.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cell.cpp b/src/cell.cpp index fc9a25e2e..35187fa9b 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -207,7 +207,7 @@ Universe::to_hdf5(hid_t universes_group) const } BoundingBox Universe::bounding_box() const { - BoundingBox bbox; + BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY}; if (cells_.size() == 0) { return {}; } else { From 72e92e9ad05e8498de26910bc9c8ef1aa747e5b6 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 18:14:57 -0500 Subject: [PATCH 27/40] More updates. Adding complex cell test to unit tests. --- include/openmc/capi.h | 4 +- src/cell.cpp | 5 +- src/material.cpp | 5 +- tests/regression_tests/complex_cell/test.py | 24 ------ tests/unit_tests/test_capi.py | 3 - tests/unit_tests/test_complex_cell_capi.py | 93 +++++++++++++++++++++ 6 files changed, 99 insertions(+), 35 deletions(-) create mode 100644 tests/unit_tests/test_complex_cell_capi.py diff --git a/include/openmc/capi.h b/include/openmc/capi.h index d81e1e722..3dd9ae697 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -14,7 +14,7 @@ extern "C" { int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n); int openmc_cell_get_id(int32_t index, int32_t* id); int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T); - int openmc_cell_get_name(int32_t index, const char*& name); + int openmc_cell_get_name(int32_t index, const char** name); int openmc_cell_set_name(int32_t index, const char* name); int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices); int openmc_cell_set_id(int32_t index, int32_t id); @@ -60,7 +60,7 @@ extern "C" { int openmc_material_set_density(int32_t index, double density, const char* units); int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density); int openmc_material_set_id(int32_t index, int32_t id); - int openmc_material_get_name(int32_t index, const char*& name); + int openmc_material_get_name(int32_t index, const char** name); int openmc_material_set_name(int32_t index, const char* name); int openmc_material_set_volume(int32_t index, double volume); int openmc_material_filter_get_bins(int32_t index, const int32_t** bins, size_t* n); diff --git a/src/cell.cpp b/src/cell.cpp index 35187fa9b..d77e0ead7 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -1085,14 +1085,13 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) //! Get the name of a cell extern "C" int -openmc_cell_get_name(int32_t index, char** name) { +openmc_cell_get_name(int32_t index, const char** name) { if (index < 0 || index >= model::cells.size()) { strcpy(openmc_err_msg, "Index in cells array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - auto name_str = model::cells[index]->name_; - strcpy(*name, name_str.c_str()); + *name = model::cells[index]->name_.c_str(); return 0; } diff --git a/src/material.cpp b/src/material.cpp index 356c2bf21..292f71b40 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1382,14 +1382,13 @@ openmc_material_set_id(int32_t index, int32_t id) } extern "C" int -openmc_material_get_name(int32_t index, char** name) { +openmc_material_get_name(int32_t index, const char** name) { if (index < 0 || index >= model::materials.size()) { strcpy(openmc_err_msg, "Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - auto name_str = model::materials[index]->name(); - strcpy(*name, name_str.c_str()); + *name = model::materials[index]->name_.data(); return 0; } diff --git a/tests/regression_tests/complex_cell/test.py b/tests/regression_tests/complex_cell/test.py index 8074c0e05..b43ccd72f 100755 --- a/tests/regression_tests/complex_cell/test.py +++ b/tests/regression_tests/complex_cell/test.py @@ -7,27 +7,3 @@ import openmc.capi def test_complex_cell(): harness = TestHarness('statepoint.10.h5') harness.main() - -def test_complex_cell_capi(): - # initialize - openmc.capi.init([]) - - inf = sys.float_info.max - - expected_boxes = { 1 : (( -4., -4., -inf), ( 4., 4., inf)), - 2 : (( -7., -7., -inf), ( 7., 7., inf)), - 3 : ((-10., -10., -inf), (10., 10., inf)), - 4 : ((-10., -10., -inf), (10., 10., inf)) } - - for cell_id, cell in openmc.capi.cells.items(): - cell_box = cell.bounding_box - - assert tuple(cell_box[0]) == expected_boxes[cell_id][0] - assert tuple(cell_box[1]) == expected_boxes[cell_id][1] - - cell_box = openmc.capi.bounding_box("Cell", cell_id) - - assert tuple(cell_box[0]) == expected_boxes[cell_id][0] - assert tuple(cell_box[1]) == expected_boxes[cell_id][1] - - openmc.capi.finalize() diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 1c1da98a3..b02e43815 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -479,9 +479,6 @@ def test_bounding_box(capi_init): llc, urc = openmc.capi.bounding_box("Surface", 5) - print(llc) - print(urc) - assert tuple(llc) == expected_llc assert tuple(urc) == expected_urc diff --git a/tests/unit_tests/test_complex_cell_capi.py b/tests/unit_tests/test_complex_cell_capi.py new file mode 100644 index 000000000..20cda54ae --- /dev/null +++ b/tests/unit_tests/test_complex_cell_capi.py @@ -0,0 +1,93 @@ +import sys + +import numpy as np +import openmc.capi + + +def test_complex_cell(run_in_tmpdir): + + openmc.reset_auto_ids() + + model = openmc.model.Model() + + u235 = openmc.Material() + u235.set_density('g/cc', 4.5) + u235.add_nuclide("U235", 1.0) + + u238 = openmc.Material() + u238.set_density('g/cc', 4.5) + u238.add_nuclide("U238", 1.0) + + zr90 = openmc.Material() + zr90.set_density('g/cc', 2.0) + zr90.add_nuclide("Zr90", 1.0) + + n14 = openmc.Material() + n14.set_density('g/cc', 0.1) + n14.add_nuclide("N14", 1.0) + + model.materials = (u235, u238, zr90, n14) + + s1 = openmc.XPlane(x0=-10.0, boundary_type='vacuum') + s2 = openmc.XPlane(x0=-7.0) + s3 = openmc.XPlane(x0=-4.0) + s4 = openmc.XPlane(x0=4.0) + s5 = openmc.XPlane(x0=7.0) + s6 = openmc.XPlane(x0=10.0, boundary_type='vacuum') + s7 = openmc.XPlane(x0=0.0) + + s11 = openmc.YPlane(y0=-10.0, boundary_type='vacuum') + s12 = openmc.YPlane(y0=-7.0) + s13 = openmc.YPlane(y0=-4.0) + s14 = openmc.YPlane(y0=4.0) + s15 = openmc.YPlane(y0=7.0) + s16 = openmc.YPlane(y0=10.0, boundary_type='vacuum') + s17 = openmc.YPlane(y0=0.0) + + c1 = openmc.Cell(fill=u235) + c1.region = +s3 & -s4 & +s13 & -s14 + + c2 = openmc.Cell(fill=u238) + c2.region = +s2 & -s5 & +s12 & -s15 & ~(+s3 & -s4 & +s13 & -s14) + + c3 = openmc.Cell(fill=zr90) + c3.region = ((+s1 & -s7 & +s17 & -s16) | (+s7 & -s6 & +s11 & -s17)) & (-s2 | +s5 | -s12 | +s15) + + c4 = openmc.Cell(fill=n14) + c4.region = ((+s1 & -s7 & +s11 & -s17) | (+s7 & -s6 & +s17 & -s16)) & ~(+s2 & -s5 & +s12 & -s15) + + model.geometry.root_universe = openmc.Universe() + model.geometry.root_universe.add_cells([c1, c2, c3, c4]) + + model.settings.batches = 10 + model.settings.inactive = 5 + model.settings.particles = 100 + model.settings.source = openmc.Source(space=openmc.stats.Box( + [-10., -10., -1.], [10., 10., 1.])) + + model.settings.verbosity = 1 + + model.export_to_xml() + + openmc.capi.finalize() + openmc.capi.init() + + inf = sys.float_info.max + + expected_boxes = { 1 : (( -4., -4., -inf), ( 4., 4., inf)), + 2 : (( -7., -7., -inf), ( 7., 7., inf)), + 3 : ((-10., -10., -inf), (10., 10., inf)), + 4 : ((-10., -10., -inf), (10., 10., inf)) } + + for cell_id, cell in openmc.capi.cells.items(): + cell_box = cell.bounding_box + + assert tuple(cell_box[0]) == expected_boxes[cell_id][0] + assert tuple(cell_box[1]) == expected_boxes[cell_id][1] + + cell_box = openmc.capi.bounding_box("Cell", cell_id) + + assert tuple(cell_box[0]) == expected_boxes[cell_id][0] + assert tuple(cell_box[1]) == expected_boxes[cell_id][1] + + openmc.capi.finalize() From 933a5c47704c5d2074e156a829f7fb00d41076d0 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 22:21:55 -0500 Subject: [PATCH 28/40] Updates based on suggestions from @paulromano. --- include/openmc/cell.h | 7 ++++--- include/openmc/surface.h | 4 ++-- openmc/capi/cell.py | 8 ++++---- openmc/capi/core.py | 10 ++++++---- openmc/capi/material.py | 16 +++++++-------- src/cell.cpp | 22 ++++++++++----------- src/geometry.cpp | 1 + tests/regression_tests/dagmc/legacy/test.py | 3 +-- 8 files changed, 36 insertions(+), 35 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 3a0c81646..752945cee 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -211,7 +211,7 @@ 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; BoundingBox bounding_box_simple() const; - BoundingBox bounding_box_complex(std::vector rpn) const; + static BoundingBox bounding_box_complex(std::vector rpn); }; //============================================================================== @@ -220,9 +220,7 @@ protected: class DAGCell : public Cell { public: - moab::DagMC* dagmc_ptr_; DAGCell(); - int32_t dag_index_; bool contains(Position r, Direction u, int32_t on_surface) const; @@ -232,6 +230,9 @@ public: BoundingBox bounding_box() const; void to_hdf5(hid_t group_id) const; + + moab::DagMC* dagmc_ptr_; //!< Pointer to DagMC instance + int32_t dag_index_; //!< DagMC index of cell }; #endif diff --git a/include/openmc/surface.h b/include/openmc/surface.h index a7d107ee3..5feb35701 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -171,8 +171,8 @@ public: void to_hdf5(hid_t group_id) const; - moab::DagMC* dagmc_ptr_; - int32_t dag_index_; + moab::DagMC* dagmc_ptr_; //!< Pointer to the DagMC instance + int32_t dag_index_; //!< DagMC index of surface }; #endif //============================================================================== diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 1388d9732..af948c92d 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -1,5 +1,5 @@ from collections.abc import Mapping, Iterable -from ctypes import byref, c_int, c_int32, c_double, c_char_p, POINTER +from ctypes import c_int, c_int32, c_double, c_char_p, POINTER from weakref import WeakValueDictionary import numpy as np @@ -111,7 +111,7 @@ class Cell(_FortranObjectWithID): @property def name(self): name = c_char_p() - _dll.openmc_cell_get_name(self._index, byref(name)) + _dll.openmc_cell_get_name(self._index, name) return name.value.decode() @name.setter @@ -185,8 +185,8 @@ class Cell(_FortranObjectWithID): @property def bounding_box(self): - llc = np.zeros((3,), dtype=float) - urc = np.zeros((3,), dtype=float) + llc = np.zeros(3) + urc = np.zeros(3) _dll.openmc_bounding_box(b'Cell', self.id, llc.ctypes.data_as(POINTER(c_double)), urc.ctypes.data_as(POINTER(c_double))) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 686dd7862..0e226d3a1 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -83,15 +83,17 @@ _dll.openmc_global_bounding_box.argtypes = [POINTER(c_double), _dll.openmc_global_bounding_box.restype = c_int _dll.openmc_global_bounding_box.errcheck = _error_handler + def global_bounding_box(): """Calculate a global bounding box for the model""" - llc = np.zeros((3,), dtype=float) - urc = np.zeros((3,), dtype=float) + llc = np.zeros(3) + urc = np.zeros(3) _dll.openmc_global_bounding_box(llc.ctypes.data_as(POINTER(c_double)), urc.ctypes.data_as(POINTER(c_double))) return llc, urc + def bounding_box(geom_type, geom_id): """Get a bounding box for a geometric object @@ -103,8 +105,8 @@ def bounding_box(geom_type, geom_id): Id of the object. Can be positive or negative for surfaces. """ geomt = c_char_p(geom_type.encode()) - llc = np.zeros((3,), dtype=float) - urc = np.zeros((3,), dtype=float) + llc = np.zeros(3) + urc = np.zeros(3) _dll.openmc_bounding_box(geomt, geom_id, llc.ctypes.data_as(POINTER(c_double)), diff --git a/openmc/capi/material.py b/openmc/capi/material.py index 3f98e9f67..f0ecac761 100644 --- a/openmc/capi/material.py +++ b/openmc/capi/material.py @@ -1,5 +1,5 @@ from collections.abc import Mapping -from ctypes import byref, c_int, c_int32, c_double, c_char_p, POINTER, c_size_t +from ctypes import c_int, c_int32, c_double, c_char_p, POINTER, c_size_t from weakref import WeakValueDictionary import numpy as np @@ -48,12 +48,12 @@ _dll.openmc_material_set_densities.errcheck = _error_handler _dll.openmc_material_set_id.argtypes = [c_int32, c_int32] _dll.openmc_material_set_id.restype = c_int _dll.openmc_material_set_id.errcheck = _error_handler -_dll.openmc_cell_get_name.argtypes = [c_int32, POINTER(c_char_p)] -_dll.openmc_cell_get_name.restype = c_int -_dll.openmc_cell_get_name.errcheck = _error_handler -_dll.openmc_cell_set_name.argtypes = [c_int32, c_char_p] -_dll.openmc_cell_set_name.restype = c_int -_dll.openmc_cell_set_name.errcheck = _error_handler +_dll.openmc_material_get_name.argtypes = [c_int32, POINTER(c_char_p)] +_dll.openmc_material_get_name.restype = c_int +_dll.openmc_material_get_name.errcheck = _error_handler +_dll.openmc_material_set_name.argtypes = [c_int32, c_char_p] +_dll.openmc_material_set_name.restype = c_int +_dll.openmc_material_set_name.errcheck = _error_handler _dll.openmc_material_set_volume.argtypes = [c_int32, c_double] _dll.openmc_material_set_volume.restype = c_int _dll.openmc_material_set_volume.errcheck = _error_handler @@ -133,7 +133,7 @@ class Material(_FortranObjectWithID): @property def name(self): name = c_char_p() - _dll.openmc_material_get_name(self._index, byref(name)) + _dll.openmc_material_get_name(self._index, name) return name.value.decode() @name.setter diff --git a/src/cell.cpp b/src/cell.cpp index d77e0ead7..35e295f78 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -1,11 +1,12 @@ #include "openmc/cell.h" +#include #include #include #include #include -#include +#include #include "openmc/capi.h" #include "openmc/constants.h" @@ -586,7 +587,7 @@ BoundingBox CSGCell::bounding_box_simple() const { return bbox; } -BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { +BoundingBox CSGCell::bounding_box_complex(std::vector rpn) { std::reverse(rpn.begin(), rpn.end()); @@ -599,7 +600,7 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { int32_t two = rpn.back(); rpn.pop_back(); // the first token should always be a surface - assert(one < OP_UNION); + Expects(one < OP_UNION); if (two >= OP_UNION) { if (two == OP_UNION) { @@ -608,14 +609,16 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { current &= model::surfaces[abs(one)-1]->bounding_box(one > 0); } } else { - // two surfaces in a row, create sub-rpn for region in parenthesis + // two surfaces in a row (left parenthesis), + // create sub-rpn for region in parenthesis std::vector subrpn; subrpn.push_back(one); subrpn.push_back(two); int32_t sone = one; int32_t stwo = two; // add until last two tokens in the sub-rpn are operators - while ((subrpn.back() < OP_UNION) || (*(subrpn.rbegin() + 1) < OP_UNION)) { + // (indicates a right parenthesis) + while (!((subrpn.back() >= OP_UNION) && (*(subrpn.rbegin() + 1) >= OP_UNION))) { subrpn.push_back(rpn.back()); rpn.pop_back(); } @@ -649,12 +652,7 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) const { } BoundingBox CSGCell::bounding_box() const { - if (simple_) { - return bounding_box_simple(); - } - else { - return bounding_box_complex(rpn_); - } + return simple_ ? bounding_box_simple() : bounding_box_complex(rpn_); } //============================================================================== @@ -1091,7 +1089,7 @@ openmc_cell_get_name(int32_t index, const char** name) { return OPENMC_E_OUT_OF_BOUNDS; } - *name = model::cells[index]->name_.c_str(); + *name = model::cells[index]->name_.data(); return 0; } diff --git a/src/geometry.cpp b/src/geometry.cpp index 43620644f..b699f2268 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -501,6 +501,7 @@ openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double } else { std::stringstream msg; msg << "Geometry type: " << gtype << " is invalid."; + set_errmsg(msg); return OPENMC_E_GEOMETRY; } diff --git a/tests/regression_tests/dagmc/legacy/test.py b/tests/regression_tests/dagmc/legacy/test.py index b6f2f55e2..d48ae9871 100644 --- a/tests/regression_tests/dagmc/legacy/test.py +++ b/tests/regression_tests/dagmc/legacy/test.py @@ -45,5 +45,4 @@ def test_dagmc(): mats = openmc.Materials([u235, water]) model.materials = mats - harness = PyAPITestHarness('statepoint.5.h5', model=model) - harness.main() + model.export_to_xml() From 5b9debe542d2bd36be5ad6630aff51bedf9da1b7 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 22:23:03 -0500 Subject: [PATCH 29/40] Update src/geometry.cpp Co-Authored-By: Paul Romano --- src/geometry.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geometry.cpp b/src/geometry.cpp index b699f2268..cc67741fd 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -519,7 +519,7 @@ openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double } extern "C" int openmc_global_bounding_box(double* llc, double* urc) { - auto bbox = model::universes[model::root_universe]->bounding_box(); + auto bbox = model::universes.at(model::root_universe)->bounding_box(); // set lower left corner values llc[0] = bbox.xmin; From 6850d4ce62a15680944bb8fe1d62dd473fd5327f Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 22:30:35 -0500 Subject: [PATCH 30/40] Apply suggestions from @paulromano Co-Authored-By: Paul Romano --- openmc/capi/core.py | 2 +- src/geometry.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 0e226d3a1..c3f5ea1b6 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -102,7 +102,7 @@ def bounding_box(geom_type, geom_id): geom_type : str Type of geometry object. One of ('surface', 'cell', 'universe') geom_id : int - Id of the object. Can be positive or negative for surfaces. + ID of the object. Can be positive or negative for surfaces. """ geomt = c_char_p(geom_type.encode()) llc = np.zeros(3) diff --git a/src/geometry.cpp b/src/geometry.cpp index cc67741fd..15118cd04 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -487,16 +487,16 @@ openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double if (gtype == "universe") { // negative ids only apply to surfaces if (id <= 0) { return OPENMC_E_GEOMETRY; } - const auto& u = model::universes[model::universe_map[id]]; + const auto& u = model::universes[model::universe_map.at(id)]; bbox = u->bounding_box(); } else if (gtype == "cell") { // negative ids only apply to surfaces if (id <= 0) { return OPENMC_E_GEOMETRY; } - const auto& c = model::cells[model::cell_map[id]]; + const auto& c = model::cells[model::cell_map.at(id)]; bbox = c->bounding_box(); } else if (gtype == "surface") { if (id == 0) { return OPENMC_E_GEOMETRY; } - const auto& s = model::surfaces[model::surface_map[abs(id)]]; + const auto& s = model::surfaces[model::surface_map.at(abs(id))]; bbox = s->bounding_box(id > 0); } else { std::stringstream msg; From c1ca3028cbfe24b47861098a81b172b205f1f9b3 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jul 2019 22:43:21 -0500 Subject: [PATCH 31/40] A few more small updates. --- include/openmc/surface.h | 2 +- src/cell.cpp | 4 ++-- src/material.cpp | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 5feb35701..d61ae07ce 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -171,7 +171,7 @@ public: void to_hdf5(hid_t group_id) const; - moab::DagMC* dagmc_ptr_; //!< Pointer to the DagMC instance + moab::DagMC* dagmc_ptr_; //!< Pointer to DagMC instance int32_t dag_index_; //!< DagMC index of surface }; #endif diff --git a/src/cell.cpp b/src/cell.cpp index 35e295f78..c3877d03b 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -1085,7 +1085,7 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) extern "C" int openmc_cell_get_name(int32_t index, const char** name) { if (index < 0 || index >= model::cells.size()) { - strcpy(openmc_err_msg, "Index in cells array is out of bounds."); + set_errmsg("Index in cells array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } @@ -1098,7 +1098,7 @@ openmc_cell_get_name(int32_t index, const char** name) { extern "C" int openmc_cell_set_name(int32_t index, const char* name) { if (index < 0 || index >= model::cells.size()) { - strcpy(openmc_err_msg, "Index in cells array is out of bounds."); + set_errmsg("Index in cells array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } diff --git a/src/material.cpp b/src/material.cpp index 292f71b40..4f3177907 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1384,7 +1384,7 @@ openmc_material_set_id(int32_t index, int32_t id) extern "C" int openmc_material_get_name(int32_t index, const char** name) { if (index < 0 || index >= model::materials.size()) { - strcpy(openmc_err_msg, "Index in materials array is out of bounds."); + set_errmsg("Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } @@ -1396,7 +1396,7 @@ openmc_material_get_name(int32_t index, const char** name) { extern "C" int openmc_material_set_name(int32_t index, const char* name) { if (index < 0 || index >= model::materials.size()) { - strcpy(openmc_err_msg, "Index in materials array is out of bounds."); + set_errmsg("Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } From dda5393c78a73126c166216686648c9a113243ca Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 19 Jul 2019 08:04:16 -0500 Subject: [PATCH 32/40] Improving error returns in bounding_box function. Appying parametrized tests. --- src/geometry.cpp | 10 ++--- tests/unit_tests/test_complex_cell_capi.py | 45 ++++++++++++---------- 2 files changed, 29 insertions(+), 26 deletions(-) diff --git a/src/geometry.cpp b/src/geometry.cpp index 15118cd04..ea90ba257 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -477,7 +477,7 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance) } extern "C" int -openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double* urc) { +openmc_bounding_box(const char* geom_type, const int32_t index, double* llc, double* urc) { BoundingBox bbox; @@ -486,23 +486,23 @@ openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double if (gtype == "universe") { // negative ids only apply to surfaces - if (id <= 0) { return OPENMC_E_GEOMETRY; } + if (id <= 0) { return OPENMC_E_INVALID_ID; } const auto& u = model::universes[model::universe_map.at(id)]; bbox = u->bounding_box(); } else if (gtype == "cell") { // negative ids only apply to surfaces - if (id <= 0) { return OPENMC_E_GEOMETRY; } + if (id <= 0) { return OPENMC_E_INVALID_ID; } const auto& c = model::cells[model::cell_map.at(id)]; bbox = c->bounding_box(); } else if (gtype == "surface") { - if (id == 0) { return OPENMC_E_GEOMETRY; } + if (id == 0) { return OPENMC_E_INVALID_IDY; } const auto& s = model::surfaces[model::surface_map.at(abs(id))]; bbox = s->bounding_box(id > 0); } else { std::stringstream msg; msg << "Geometry type: " << gtype << " is invalid."; set_errmsg(msg); - return OPENMC_E_GEOMETRY; + return OPENMC_E_INVALID_TYPE; } // set lower left corner values diff --git a/tests/unit_tests/test_complex_cell_capi.py b/tests/unit_tests/test_complex_cell_capi.py index 20cda54ae..ebc345f8f 100644 --- a/tests/unit_tests/test_complex_cell_capi.py +++ b/tests/unit_tests/test_complex_cell_capi.py @@ -2,9 +2,10 @@ import sys import numpy as np import openmc.capi +import pytest - -def test_complex_cell(run_in_tmpdir): +@pytest.fixture(autouse=True) +def complex_cell(run_in_tmpdir): openmc.reset_auto_ids() @@ -51,10 +52,12 @@ def test_complex_cell(run_in_tmpdir): c2.region = +s2 & -s5 & +s12 & -s15 & ~(+s3 & -s4 & +s13 & -s14) c3 = openmc.Cell(fill=zr90) - c3.region = ((+s1 & -s7 & +s17 & -s16) | (+s7 & -s6 & +s11 & -s17)) & (-s2 | +s5 | -s12 | +s15) + c3.region = ((+s1 & -s7 & +s17 & -s16) | (+s7 & -s6 & +s11 & -s17)) \ + & (-s2 | +s5 | -s12 | +s15) c4 = openmc.Cell(fill=n14) - c4.region = ((+s1 & -s7 & +s11 & -s17) | (+s7 & -s6 & +s17 & -s16)) & ~(+s2 & -s5 & +s12 & -s15) + c4.region = ((+s1 & -s7 & +s11 & -s17) | (+s7 & -s6 & +s17 & -s16)) & \ + ~(+s2 & -s5 & +s12 & -s15) model.geometry.root_universe = openmc.Universe() model.geometry.root_universe.add_cells([c1, c2, c3, c4]) @@ -72,22 +75,22 @@ def test_complex_cell(run_in_tmpdir): openmc.capi.finalize() openmc.capi.init() - inf = sys.float_info.max - - expected_boxes = { 1 : (( -4., -4., -inf), ( 4., 4., inf)), - 2 : (( -7., -7., -inf), ( 7., 7., inf)), - 3 : ((-10., -10., -inf), (10., 10., inf)), - 4 : ((-10., -10., -inf), (10., 10., inf)) } - - for cell_id, cell in openmc.capi.cells.items(): - cell_box = cell.bounding_box - - assert tuple(cell_box[0]) == expected_boxes[cell_id][0] - assert tuple(cell_box[1]) == expected_boxes[cell_id][1] - - cell_box = openmc.capi.bounding_box("Cell", cell_id) - - assert tuple(cell_box[0]) == expected_boxes[cell_id][0] - assert tuple(cell_box[1]) == expected_boxes[cell_id][1] + yield openmc.capi.finalize() + +inf = sys.float_info.max + +expected_results = ( (1, (( -4., -4., -inf), ( 4., 4., inf))), + (2, (( -7., -7., -inf), ( 7., 7., inf))), + (3, ((-10., -10., -inf), (10., 10., inf))), + (4, ((-10., -10., -inf), (10., 10., inf))) ) +@pytest.mark.parametrize("cell_id,expected_box", expected_results) +def test_cell_box(cell_id, expected_box): + cell_box = openmc.capi.cells[cell_id].bounding_box + assert tuple(cell_box[0]) == expected_box[0] + assert tuple(cell_box[1]) == expected_box[1] + + cell_box = openmc.capi.bounding_box("Cell", cell_id) + assert tuple(cell_box[0]) == expected_box[0] + assert tuple(cell_box[1]) == expected_box[1] From 694bd0cf15efd3467827dbf5247d6afd66a0244a Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 19 Jul 2019 08:53:52 -0500 Subject: [PATCH 33/40] Chaning back to id in bounding_box for now. --- src/geometry.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/geometry.cpp b/src/geometry.cpp index ea90ba257..91e2fc245 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -477,7 +477,7 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance) } extern "C" int -openmc_bounding_box(const char* geom_type, const int32_t index, double* llc, double* urc) { +openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double* urc) { BoundingBox bbox; @@ -495,7 +495,7 @@ openmc_bounding_box(const char* geom_type, const int32_t index, double* llc, dou const auto& c = model::cells[model::cell_map.at(id)]; bbox = c->bounding_box(); } else if (gtype == "surface") { - if (id == 0) { return OPENMC_E_INVALID_IDY; } + if (id == 0) { return OPENMC_E_INVALID_ID; } const auto& s = model::surfaces[model::surface_map.at(abs(id))]; bbox = s->bounding_box(id > 0); } else { From f2199af0eba57a001d196d96bac0ad66d1d1d842 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 19 Jul 2019 09:46:47 -0500 Subject: [PATCH 34/40] Cleaning up error checking in bounding_box. Updating expected exceptions in the capi tests. --- src/geometry.cpp | 19 ++++++++++++++----- tests/unit_tests/test_capi.py | 6 +++--- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/src/geometry.cpp b/src/geometry.cpp index 91e2fc245..ae48a0acb 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -484,18 +484,27 @@ openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double std::string gtype(geom_type); to_lower(gtype); + // negative ids only for surfaces + if (id < 0 && gtype != "surface") { + std::stringstream msg; + msg << "Negative ID " << id << " passed to bounding box for " + << "non-surface geom type '" << geom_type << "'" << std::endl; + set_errmsg(msg); + return OPENMC_E_INVALID_ID; + } + // id should never be zero + if (id == 0) { + set_errmsg("Invalid ID 0 for surface bounding box."); + return OPENMC_E_INVALID_ID; + } + if (gtype == "universe") { - // negative ids only apply to surfaces - if (id <= 0) { return OPENMC_E_INVALID_ID; } const auto& u = model::universes[model::universe_map.at(id)]; bbox = u->bounding_box(); } else if (gtype == "cell") { - // negative ids only apply to surfaces - if (id <= 0) { return OPENMC_E_INVALID_ID; } const auto& c = model::cells[model::cell_map.at(id)]; bbox = c->bounding_box(); } else if (gtype == "surface") { - if (id == 0) { return OPENMC_E_INVALID_ID; } const auto& s = model::surfaces[model::surface_map.at(abs(id))]; bbox = s->bounding_box(id > 0); } else { diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index b02e43815..aa140a121 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -508,13 +508,13 @@ def test_bounding_box(capi_init): assert tuple(urc) == expected_urc # make sure that proper assertions are raised - with pytest.raises(openmc.exceptions.GeometryError): + with pytest.raises(openmc.exceptions.InvalidIDError): openmc.capi.bounding_box("Cell", -1) - with pytest.raises(openmc.exceptions.GeometryError): + with pytest.raises(openmc.exceptions.InvalidIDError): openmc.capi.bounding_box("Surface", 0) - with pytest.raises(openmc.exceptions.GeometryError): + with pytest.raises(openmc.exceptions.InvalidTypeError): openmc.capi.bounding_box("Region", 1) From 331e198617d7f68fab1ba4af5da0143d07b4e025 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 19 Jul 2019 10:03:22 -0500 Subject: [PATCH 35/40] Adding test for name setting and relying on std::string's char* constructor. --- src/cell.cpp | 3 +-- src/material.cpp | 1 - tests/unit_tests/test_capi.py | 4 ++++ 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index c3877d03b..7cc74ed9f 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -1102,8 +1102,7 @@ openmc_cell_set_name(int32_t index, const char* name) { return OPENMC_E_OUT_OF_BOUNDS; } - std::string name_str(name); - model::cells[index]->set_name(name_str); + model::cells[index]->set_name(name); return 0; } diff --git a/src/material.cpp b/src/material.cpp index 4f3177907..699649d31 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1400,7 +1400,6 @@ openmc_material_set_name(int32_t index, const char* name) { return OPENMC_E_OUT_OF_BOUNDS; } - std::string name_str(name); model::materials[index]->set_name(name); return 0; diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index aa140a121..3fdd4db66 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -75,6 +75,8 @@ def test_cell(capi_init): cell.fill = openmc.capi.materials[1] assert str(cell) == 'Cell[0]' assert cell.name == "Fuel" + cell.name = "Not fuel" + assert cell.name == "Not fuel" def test_cell_temperature(capi_init): cell = openmc.capi.cells[1] @@ -124,6 +126,8 @@ def test_material(capi_init): m.set_density(0.1, 'g/cm3') assert m.density == pytest.approx(0.1) assert m.name == "Hot borated water" + m.name = "Not hot borated water" + assert m.name == "Not hot borated water" def test_material_add_nuclide(capi_init): m = openmc.capi.materials[3] From d820e6777aef9556c355e183be64032ad269f5b8 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 23 Jul 2019 14:33:35 -0500 Subject: [PATCH 36/40] Apply suggestions from code review Accessing cell/material names through accessor. Returning a const reference from `name()` accessor in both cases. Co-Authored-By: Paul Romano --- include/openmc/cell.h | 2 +- include/openmc/material.h | 2 +- src/cell.cpp | 2 +- src/material.cpp | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 752945cee..602da98e7 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -136,7 +136,7 @@ public: //! Get the name of a cell //! \return Cell name - std::string name() const { return name_; }; + const std::string& name() const { return name_; }; //! Set the temperature of a cell instance //! \param[in] name Cell name diff --git a/include/openmc/material.h b/include/openmc/material.h index aa0a4bad0..652f3db8e 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -94,7 +94,7 @@ public: //! Get name //! \return Material name - std::string name() const { return name_; } + const std::string& name() const { return name_; } //! Set name void set_name(const std::string& name) { name_ = name; } diff --git a/src/cell.cpp b/src/cell.cpp index 7cc74ed9f..49d8cb1b4 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -1089,7 +1089,7 @@ openmc_cell_get_name(int32_t index, const char** name) { return OPENMC_E_OUT_OF_BOUNDS; } - *name = model::cells[index]->name_.data(); + *name = model::cells[index]->name().data(); return 0; } diff --git a/src/material.cpp b/src/material.cpp index 699649d31..685fd1da2 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1388,7 +1388,7 @@ openmc_material_get_name(int32_t index, const char** name) { return OPENMC_E_OUT_OF_BOUNDS; } - *name = model::materials[index]->name_.data(); + *name = model::materials[index]->name().data(); return 0; } From 3439797e322a78f88bce137ad77835e3c2fb3706 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 24 Jul 2019 03:02:10 -0500 Subject: [PATCH 37/40] Update for case where the entire region is a complement. Addition of a cell and redefinition of another to test more robustly. --- include/openmc/cell.h | 1 + src/cell.cpp | 23 +++++++++++++++++----- tests/unit_tests/test_complex_cell_capi.py | 15 +++++++------- 3 files changed, 27 insertions(+), 12 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 602da98e7..db8d60ed4 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -212,6 +212,7 @@ protected: bool contains_complex(Position r, Direction u, int32_t on_surface) const; BoundingBox bounding_box_simple() const; static BoundingBox bounding_box_complex(std::vector rpn); + static void apply_demorgan(std::vector& rpn); }; //============================================================================== diff --git a/src/cell.cpp b/src/cell.cpp index 49d8cb1b4..510002fa9 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -587,8 +587,25 @@ BoundingBox CSGCell::bounding_box_simple() const { return bbox; } +void CSGCell::apply_demorgan(std::vector& rpn) { + for (auto& token : rpn) { + if (token < OP_UNION) { token *= -1; } + else if (token == OP_UNION) { token = OP_INTERSECTION; } + else if (token == OP_INTERSECTION) { token = OP_UNION; } + } +} + BoundingBox CSGCell::bounding_box_complex(std::vector rpn) { + // if the last operator is a complement op, there is no + // sub-region that the complement connects to. This indicates + // that the entire region is a complement and we can apply + // De Morgan's laws immediately + if ((rpn.back() == OP_COMPLEMENT)) { + rpn.pop_back(); + apply_demorgan(rpn); + } + std::reverse(rpn.begin(), rpn.end()); BoundingBox current = model::surfaces[abs(rpn.back()) - 1]->bounding_box(rpn.back() > 0); @@ -626,11 +643,7 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) { // handle complement case using De Morgan's laws if (subrpn.back() == OP_COMPLEMENT) { subrpn.pop_back(); - for (auto& token : subrpn) { - if (token < OP_UNION) { token *= -1; } - else if (token == OP_UNION) { token = OP_INTERSECTION; } - else if (token == OP_INTERSECTION) { token = OP_UNION; } - } + apply_demorgan(subrpn); subrpn.push_back(rpn.back()); rpn.pop_back(); } diff --git a/tests/unit_tests/test_complex_cell_capi.py b/tests/unit_tests/test_complex_cell_capi.py index ebc345f8f..451233d64 100644 --- a/tests/unit_tests/test_complex_cell_capi.py +++ b/tests/unit_tests/test_complex_cell_capi.py @@ -46,7 +46,7 @@ def complex_cell(run_in_tmpdir): s17 = openmc.YPlane(y0=0.0) c1 = openmc.Cell(fill=u235) - c1.region = +s3 & -s4 & +s13 & -s14 + c1.region = ~(-s3 | +s4 | ~(+s13 & -s14)) c2 = openmc.Cell(fill=u238) c2.region = +s2 & -s5 & +s12 & -s15 & ~(+s3 & -s4 & +s13 & -s14) @@ -59,8 +59,11 @@ def complex_cell(run_in_tmpdir): c4.region = ((+s1 & -s7 & +s11 & -s17) | (+s7 & -s6 & +s17 & -s16)) & \ ~(+s2 & -s5 & +s12 & -s15) + c5 = openmc.Cell(fill=n14) + c5.region = ~(+s1 & -s6 & +s11 & -s16) + model.geometry.root_universe = openmc.Universe() - model.geometry.root_universe.add_cells([c1, c2, c3, c4]) + model.geometry.root_universe.add_cells([c1, c2, c3, c4, c5]) model.settings.batches = 10 model.settings.inactive = 5 @@ -84,13 +87,11 @@ inf = sys.float_info.max expected_results = ( (1, (( -4., -4., -inf), ( 4., 4., inf))), (2, (( -7., -7., -inf), ( 7., 7., inf))), (3, ((-10., -10., -inf), (10., 10., inf))), - (4, ((-10., -10., -inf), (10., 10., inf))) ) + (4, ((-10., -10., -inf), (10., 10., inf))), + (5, ((-inf, -inf, -inf), (inf, inf, inf))) ) @pytest.mark.parametrize("cell_id,expected_box", expected_results) def test_cell_box(cell_id, expected_box): + print("Cell {}".format(cell_id)) cell_box = openmc.capi.cells[cell_id].bounding_box assert tuple(cell_box[0]) == expected_box[0] assert tuple(cell_box[1]) == expected_box[1] - - cell_box = openmc.capi.bounding_box("Cell", cell_id) - assert tuple(cell_box[0]) == expected_box[0] - assert tuple(cell_box[1]) == expected_box[1] From 3be6520671b7fc2fda51446edeffd5aca0f13225 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 24 Jul 2019 12:49:02 -0500 Subject: [PATCH 38/40] Removing print statement. --- tests/unit_tests/test_complex_cell_capi.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/unit_tests/test_complex_cell_capi.py b/tests/unit_tests/test_complex_cell_capi.py index 451233d64..e5983d2d5 100644 --- a/tests/unit_tests/test_complex_cell_capi.py +++ b/tests/unit_tests/test_complex_cell_capi.py @@ -91,7 +91,6 @@ expected_results = ( (1, (( -4., -4., -inf), ( 4., 4., inf))), (5, ((-inf, -inf, -inf), (inf, inf, inf))) ) @pytest.mark.parametrize("cell_id,expected_box", expected_results) def test_cell_box(cell_id, expected_box): - print("Cell {}".format(cell_id)) cell_box = openmc.capi.cells[cell_id].bounding_box assert tuple(cell_box[0]) == expected_box[0] assert tuple(cell_box[1]) == expected_box[1] From 5ad4f94f3adcb2653b7c91d8c6c636e4238d1e96 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 24 Jul 2019 13:10:03 -0500 Subject: [PATCH 39/40] Removing CAPI function for bounding box of any type. Moving to a cell-only function. Tests are adjusted as needed. --- include/openmc/capi.h | 2 +- openmc/capi/cell.py | 7 ++++- openmc/capi/core.py | 27 ------------------- src/cell.cpp | 28 ++++++++++++++++--- src/geometry.cpp | 51 ----------------------------------- tests/unit_tests/test_capi.py | 46 ------------------------------- 6 files changed, 32 insertions(+), 129 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 3dd9ae697..ded41f738 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -31,7 +31,7 @@ extern "C" { int openmc_filter_set_id(int32_t index, int32_t id); int openmc_finalize(); int openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance); - int openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double* urc); + int openmc_cell_bounding_box(const int32_t index, double* llc, double* urc); int openmc_global_bounding_box(double* llc, double* urc); int openmc_fission_bank(void** ptr, int64_t* n); int openmc_get_cell_index(int32_t id, int32_t* index); diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index af948c92d..de8ee9f43 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -49,6 +49,11 @@ _dll.openmc_get_cell_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_cell_index.restype = c_int _dll.openmc_get_cell_index.errcheck = _error_handler _dll.cells_size.restype = c_int +_dll.openmc_cell_bounding_box.argtypes = [c_int, + POINTER(c_double), + POINTER(c_double)] +_dll.openmc_cell_bounding_box.restype = c_int +_dll.openmc_cell_bounding_box.errcheck = _error_handler class Cell(_FortranObjectWithID): @@ -187,7 +192,7 @@ class Cell(_FortranObjectWithID): def bounding_box(self): llc = np.zeros(3) urc = np.zeros(3) - _dll.openmc_bounding_box(b'Cell', self.id, + _dll.openmc_cell_bounding_box(self._index, llc.ctypes.data_as(POINTER(c_double)), urc.ctypes.data_as(POINTER(c_double))) return llc, urc diff --git a/openmc/capi/core.py b/openmc/capi/core.py index c3f5ea1b6..cd7401a23 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -73,11 +73,6 @@ _dll.openmc_simulation_finalize.errcheck = _error_handler _dll.openmc_statepoint_write.argtypes = [c_char_p, POINTER(c_bool)] _dll.openmc_statepoint_write.restype = c_int _dll.openmc_statepoint_write.errcheck = _error_handler -_dll.openmc_bounding_box.argtypes = [c_char_p, c_int, - POINTER(c_double), - POINTER(c_double)] -_dll.openmc_bounding_box.restype = c_int -_dll.openmc_bounding_box.errcheck = _error_handler _dll.openmc_global_bounding_box.argtypes = [POINTER(c_double), POINTER(c_double)] _dll.openmc_global_bounding_box.restype = c_int @@ -93,28 +88,6 @@ def global_bounding_box(): return llc, urc - -def bounding_box(geom_type, geom_id): - """Get a bounding box for a geometric object - - Parameters - ---------- - geom_type : str - Type of geometry object. One of ('surface', 'cell', 'universe') - geom_id : int - ID of the object. Can be positive or negative for surfaces. - """ - geomt = c_char_p(geom_type.encode()) - llc = np.zeros(3) - urc = np.zeros(3) - _dll.openmc_bounding_box(geomt, - geom_id, - llc.ctypes.data_as(POINTER(c_double)), - urc.ctypes.data_as(POINTER(c_double))) - - return llc, urc - - def calculate_volumes(): """Run stochastic volume calculation""" _dll.openmc_calculate_volumes() diff --git a/src/cell.cpp b/src/cell.cpp index 510002fa9..3c079f7ae 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -587,6 +587,7 @@ BoundingBox CSGCell::bounding_box_simple() const { return bbox; } + void CSGCell::apply_demorgan(std::vector& rpn) { for (auto& token : rpn) { if (token < OP_UNION) { token *= -1; } @@ -601,11 +602,12 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) { // sub-region that the complement connects to. This indicates // that the entire region is a complement and we can apply // De Morgan's laws immediately - if ((rpn.back() == OP_COMPLEMENT)) { + if (rpn.back() == OP_COMPLEMENT) { rpn.pop_back(); apply_demorgan(rpn); } + // reverse the rpn to make popping easier std::reverse(rpn.begin(), rpn.end()); BoundingBox current = model::surfaces[abs(rpn.back()) - 1]->bounding_box(rpn.back() > 0); @@ -631,8 +633,6 @@ BoundingBox CSGCell::bounding_box_complex(std::vector rpn) { std::vector subrpn; subrpn.push_back(one); subrpn.push_back(two); - int32_t sone = one; - int32_t stwo = two; // add until last two tokens in the sub-rpn are operators // (indicates a right parenthesis) while (!((subrpn.back() >= OP_UNION) && (*(subrpn.rbegin() + 1) >= OP_UNION))) { @@ -1094,6 +1094,28 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) return 0; } +//! Get the bounding box of a cell +extern "C" int +openmc_cell_bounding_box(const int32_t index, double* llc, double* urc) { + + BoundingBox bbox; + + const auto& c = model::cells[index]; + bbox = c->bounding_box(); + + // set lower left corner values + llc[0] = bbox.xmin; + llc[1] = bbox.ymin; + llc[2] = bbox.zmin; + + // set upper right corner values + urc[0] = bbox.xmax; + urc[1] = bbox.ymax; + urc[2] = bbox.zmax; + + return 0; +} + //! Get the name of a cell extern "C" int openmc_cell_get_name(int32_t index, const char** name) { diff --git a/src/geometry.cpp b/src/geometry.cpp index ae48a0acb..556d10650 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -476,57 +476,6 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance) return 0; } -extern "C" int -openmc_bounding_box(const char* geom_type, const int32_t id, double* llc, double* urc) { - - BoundingBox bbox; - - std::string gtype(geom_type); - to_lower(gtype); - - // negative ids only for surfaces - if (id < 0 && gtype != "surface") { - std::stringstream msg; - msg << "Negative ID " << id << " passed to bounding box for " - << "non-surface geom type '" << geom_type << "'" << std::endl; - set_errmsg(msg); - return OPENMC_E_INVALID_ID; - } - // id should never be zero - if (id == 0) { - set_errmsg("Invalid ID 0 for surface bounding box."); - return OPENMC_E_INVALID_ID; - } - - if (gtype == "universe") { - const auto& u = model::universes[model::universe_map.at(id)]; - bbox = u->bounding_box(); - } else if (gtype == "cell") { - const auto& c = model::cells[model::cell_map.at(id)]; - bbox = c->bounding_box(); - } else if (gtype == "surface") { - const auto& s = model::surfaces[model::surface_map.at(abs(id))]; - bbox = s->bounding_box(id > 0); - } else { - std::stringstream msg; - msg << "Geometry type: " << gtype << " is invalid."; - set_errmsg(msg); - return OPENMC_E_INVALID_TYPE; - } - - // set lower left corner values - llc[0] = bbox.xmin; - llc[1] = bbox.ymin; - llc[2] = bbox.zmin; - - // set upper right corner values - urc[0] = bbox.xmax; - urc[1] = bbox.ymax; - urc[2] = bbox.zmax; - - return 0; -} - extern "C" int openmc_global_bounding_box(double* llc, double* urc) { auto bbox = model::universes.at(model::root_universe)->bounding_box(); diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 3fdd4db66..d1c578ba0 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -475,52 +475,6 @@ def test_position(capi_init): assert tuple(pos) == (1.3, 2.3, 3.3) -def test_bounding_box(capi_init): - inf = sys.float_info.max - - expected_llc = (-inf, -0.63, -inf) - expected_urc = (inf, inf, inf) - - llc, urc = openmc.capi.bounding_box("Surface", 5) - - assert tuple(llc) == expected_llc - assert tuple(urc) == expected_urc - - - expected_llc = (-inf, -inf, -inf) - expected_urc = (inf, -0.63, inf) - - llc, urc = openmc.capi.bounding_box("Surface", -5) - - assert tuple(llc) == expected_llc - assert tuple(urc) == expected_urc - - expected_llc = (-0.39218, -0.39218, -inf) - expected_urc = (0.39218, 0.39218, inf) - - llc, urc = openmc.capi.bounding_box("Cell", 1) - - assert tuple(llc) == expected_llc - assert tuple(urc) == expected_urc - - expected_llc = (-0.45720, -0.45720, -inf) - expected_urc = (0.45720, 0.45720, inf) - - llc, urc = openmc.capi.bounding_box("Cell", 2) - - assert tuple(llc) == expected_llc - assert tuple(urc) == expected_urc - - # make sure that proper assertions are raised - with pytest.raises(openmc.exceptions.InvalidIDError): - openmc.capi.bounding_box("Cell", -1) - - with pytest.raises(openmc.exceptions.InvalidIDError): - openmc.capi.bounding_box("Surface", 0) - - with pytest.raises(openmc.exceptions.InvalidTypeError): - openmc.capi.bounding_box("Region", 1) - def test_global_bounding_box(capi_init): inf = sys.float_info.max From d5eaa3b5b768a67a3c80b4d033ecc7dbc6695db3 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 24 Jul 2019 14:36:55 -0500 Subject: [PATCH 40/40] Converting double max to np.inf on the Python side for bounding boxes. --- openmc/capi/cell.py | 8 ++++++++ openmc/capi/core.py | 7 +++++++ tests/unit_tests/test_capi.py | 7 ++----- tests/unit_tests/test_complex_cell_capi.py | 18 ++++++++++-------- 4 files changed, 27 insertions(+), 13 deletions(-) diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index de8ee9f43..784ebd087 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -1,3 +1,5 @@ +import sys + from collections.abc import Mapping, Iterable from ctypes import c_int, c_int32, c_double, c_char_p, POINTER from weakref import WeakValueDictionary @@ -190,11 +192,17 @@ class Cell(_FortranObjectWithID): @property def bounding_box(self): + inf = sys.float_info.max llc = np.zeros(3) urc = np.zeros(3) _dll.openmc_cell_bounding_box(self._index, llc.ctypes.data_as(POINTER(c_double)), urc.ctypes.data_as(POINTER(c_double))) + llc[llc == inf] = np.inf + urc[urc == inf] = np.inf + llc[llc == -inf] = -np.inf + urc[urc == -inf] = -np.inf + return llc, urc class _CellMapping(Mapping): diff --git a/openmc/capi/core.py b/openmc/capi/core.py index cd7401a23..a470f0665 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -1,3 +1,5 @@ +import sys + from contextlib import contextmanager from ctypes import (CDLL, c_bool, c_int, c_int32, c_int64, c_double, c_char_p, c_char, POINTER, Structure, c_void_p, create_string_buffer) @@ -81,10 +83,15 @@ _dll.openmc_global_bounding_box.errcheck = _error_handler def global_bounding_box(): """Calculate a global bounding box for the model""" + inf = sys.float_info.max llc = np.zeros(3) urc = np.zeros(3) _dll.openmc_global_bounding_box(llc.ctypes.data_as(POINTER(c_double)), urc.ctypes.data_as(POINTER(c_double))) + llc[llc == inf] = np.inf + urc[urc == inf] = np.inf + llc[llc == -inf] = -np.inf + urc[urc == -inf] = -np.inf return llc, urc diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index d1c578ba0..6953186eb 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -1,6 +1,5 @@ from collections.abc import Mapping import os -import sys import numpy as np import pytest @@ -477,10 +476,8 @@ def test_position(capi_init): def test_global_bounding_box(capi_init): - inf = sys.float_info.max - - expected_llc = (-0.63, -0.63, -inf) - expected_urc = (0.63, 0.63, inf) + expected_llc = (-0.63, -0.63, -np.inf) + expected_urc = (0.63, 0.63, np.inf) llc, urc = openmc.capi.global_bounding_box() diff --git a/tests/unit_tests/test_complex_cell_capi.py b/tests/unit_tests/test_complex_cell_capi.py index e5983d2d5..365ce4808 100644 --- a/tests/unit_tests/test_complex_cell_capi.py +++ b/tests/unit_tests/test_complex_cell_capi.py @@ -1,5 +1,3 @@ -import sys - import numpy as np import openmc.capi import pytest @@ -82,13 +80,17 @@ def complex_cell(run_in_tmpdir): openmc.capi.finalize() -inf = sys.float_info.max -expected_results = ( (1, (( -4., -4., -inf), ( 4., 4., inf))), - (2, (( -7., -7., -inf), ( 7., 7., inf))), - (3, ((-10., -10., -inf), (10., 10., inf))), - (4, ((-10., -10., -inf), (10., 10., inf))), - (5, ((-inf, -inf, -inf), (inf, inf, inf))) ) +expected_results = ( (1, (( -4., -4., -np.inf), + ( 4., 4., np.inf))), + (2, (( -7., -7., -np.inf), + ( 7., 7., np.inf))), + (3, ((-10., -10., -np.inf), + ( 10., 10., np.inf))), + (4, ((-10., -10., -np.inf), + ( 10., 10., np.inf))), + (5, ((-np.inf, -np.inf, -np.inf), + ( np.inf, np.inf, np.inf))) ) @pytest.mark.parametrize("cell_id,expected_box", expected_results) def test_cell_box(cell_id, expected_box): cell_box = openmc.capi.cells[cell_id].bounding_box