mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge pull request #1295 from pshriwise/bounding_boxes
Bounding box updates
This commit is contained in:
Paul Romano
GitHub
commit
5178272ca8
15 changed files with 540 additions and 49 deletions
|
|
@ -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);
|
||||
|
|
@ -29,6 +31,8 @@ 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_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);
|
||||
int openmc_get_filter_index(int32_t id, int32_t* index);
|
||||
|
|
@ -56,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);
|
||||
|
|
|
|||
|
|
@ -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<UniversePartitioner> partitioner_;
|
||||
};
|
||||
|
||||
|
|
@ -114,6 +116,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
|
||||
|
||||
|
|
@ -129,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
|
||||
const std::string& name() const { return name_; };
|
||||
|
||||
//! Set the temperature of a cell instance
|
||||
//! \param[in] name Cell name
|
||||
void set_name(const std::string& name) { name_ = name; };
|
||||
|
||||
//----------------------------------------------------------------------------
|
||||
// Data members
|
||||
|
||||
|
|
@ -192,9 +205,14 @@ 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;
|
||||
BoundingBox bounding_box_simple() const;
|
||||
static BoundingBox bounding_box_complex(std::vector<int32_t> rpn);
|
||||
static void apply_demorgan(std::vector<int32_t>& rpn);
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -203,16 +221,19 @@ 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;
|
||||
|
||||
std::pair<double, int32_t>
|
||||
distance(Position r, Direction u, int32_t on_surface) const;
|
||||
|
||||
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
|
||||
|
||||
|
|
|
|||
|
|
@ -92,6 +92,13 @@ public:
|
|||
//! \return Density in [g/cm^3]
|
||||
double density_gpcc() const { return density_gpcc_; }
|
||||
|
||||
//! Get name
|
||||
//! \return Material name
|
||||
const 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
|
||||
|
|
|
|||
|
|
@ -43,12 +43,46 @@ 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;
|
||||
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -105,6 +139,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(bool pos_side) const { return {}; }
|
||||
};
|
||||
|
||||
class CSGSurface : public Surface
|
||||
|
|
@ -126,18 +162,17 @@ 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;
|
||||
Direction normal(Position r) const;
|
||||
Direction reflect(Position r, Direction u) const;
|
||||
//! Get the bounding box of this surface.
|
||||
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 surface
|
||||
};
|
||||
#endif
|
||||
//==============================================================================
|
||||
|
|
@ -166,8 +201,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() const = 0;
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -186,7 +219,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_;
|
||||
};
|
||||
|
|
@ -207,7 +240,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_;
|
||||
};
|
||||
|
|
@ -228,7 +261,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_;
|
||||
};
|
||||
|
|
@ -249,7 +282,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() const;
|
||||
|
||||
double A_, B_, C_, D_;
|
||||
};
|
||||
|
|
@ -269,6 +301,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(bool pos_side) const;
|
||||
|
||||
double y0_, z0_, radius_;
|
||||
};
|
||||
|
|
@ -288,6 +321,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(bool pos_side) const;
|
||||
|
||||
double x0_, z0_, radius_;
|
||||
};
|
||||
|
|
@ -307,6 +341,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(bool pos_side) const;
|
||||
|
||||
double x0_, y0_, radius_;
|
||||
};
|
||||
|
|
@ -326,6 +361,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(bool pos_side) const;
|
||||
|
||||
double x0_, y0_, z0_, radius_;
|
||||
};
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -28,6 +30,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
|
||||
|
|
@ -43,6 +51,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):
|
||||
|
|
@ -102,6 +115,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, 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()
|
||||
|
|
@ -166,6 +190,20 @@ class Cell(_FortranObjectWithID):
|
|||
|
||||
_dll.openmc_cell_set_temperature(self._index, T, instance)
|
||||
|
||||
@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):
|
||||
def __getitem__(self, key):
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
@ -73,8 +75,26 @@ _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_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():
|
||||
"""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
|
||||
|
||||
def calculate_volumes():
|
||||
"""Run stochastic volume calculation"""
|
||||
_dll.openmc_calculate_volumes()
|
||||
|
|
|
|||
|
|
@ -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_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
|
||||
|
|
@ -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, 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()
|
||||
|
|
|
|||
166
src/cell.cpp
166
src/cell.cpp
|
|
@ -1,11 +1,12 @@
|
|||
|
||||
#include "openmc/cell.h"
|
||||
|
||||
#include <cctype>
|
||||
#include <cmath>
|
||||
#include <sstream>
|
||||
#include <set>
|
||||
#include <string>
|
||||
#include <cctype>
|
||||
#include <gsl/gsl>
|
||||
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/constants.h"
|
||||
|
|
@ -118,7 +119,6 @@ generate_rpn(int32_t cell_id, std::vector<int32_t> infix)
|
|||
if (token < OP_UNION) {
|
||||
// If token is not an operator, add it to output
|
||||
rpn.push_back(token);
|
||||
|
||||
} else if (token < OP_RIGHT_PAREN) {
|
||||
// Regular operators union, intersection, complement
|
||||
while (stack.size() > 0) {
|
||||
|
|
@ -158,7 +158,6 @@ generate_rpn(int32_t cell_id, std::vector<int32_t> infix)
|
|||
<< cell_id;
|
||||
fatal_error(err_msg);
|
||||
}
|
||||
|
||||
rpn.push_back(stack.back());
|
||||
stack.pop_back();
|
||||
}
|
||||
|
|
@ -208,6 +207,19 @@ Universe::to_hdf5(hid_t universes_group) const
|
|||
close_group(group);
|
||||
}
|
||||
|
||||
BoundingBox Universe::bounding_box() const {
|
||||
BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
|
||||
if (cells_.size() == 0) {
|
||||
return {};
|
||||
} else {
|
||||
for (const auto& cell : cells_) {
|
||||
auto& c = model::cells[cell];
|
||||
bbox |= c->bounding_box();
|
||||
}
|
||||
}
|
||||
return bbox;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Cell implementation
|
||||
//==============================================================================
|
||||
|
|
@ -567,6 +579,95 @@ CSGCell::to_hdf5(hid_t cell_group) const
|
|||
close_group(group);
|
||||
}
|
||||
|
||||
BoundingBox CSGCell::bounding_box_simple() const {
|
||||
BoundingBox bbox;
|
||||
for (int32_t token : rpn_) {
|
||||
bbox &= model::surfaces[abs(token)-1]->bounding_box(token > 0);
|
||||
}
|
||||
return bbox;
|
||||
}
|
||||
|
||||
|
||||
void CSGCell::apply_demorgan(std::vector<int32_t>& 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<int32_t> 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);
|
||||
}
|
||||
|
||||
// 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);
|
||||
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
|
||||
Expects(one < OP_UNION);
|
||||
|
||||
if (two >= OP_UNION) {
|
||||
if (two == OP_UNION) {
|
||||
current |= model::surfaces[abs(one)-1]->bounding_box(one > 0);
|
||||
} else if (two == OP_INTERSECTION) {
|
||||
current &= model::surfaces[abs(one)-1]->bounding_box(one > 0);
|
||||
}
|
||||
} else {
|
||||
// two surfaces in a row (left parenthesis),
|
||||
// create sub-rpn for region in parenthesis
|
||||
std::vector<int32_t> subrpn;
|
||||
subrpn.push_back(one);
|
||||
subrpn.push_back(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))) {
|
||||
subrpn.push_back(rpn.back());
|
||||
rpn.pop_back();
|
||||
}
|
||||
|
||||
// handle complement case using De Morgan's laws
|
||||
if (subrpn.back() == OP_COMPLEMENT) {
|
||||
subrpn.pop_back();
|
||||
apply_demorgan(subrpn);
|
||||
subrpn.push_back(rpn.back());
|
||||
rpn.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 |= sub_box;
|
||||
} else if (op == OP_INTERSECTION) {
|
||||
current &= sub_box;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return current;
|
||||
}
|
||||
|
||||
BoundingBox CSGCell::bounding_box() const {
|
||||
return simple_ ? bounding_box_simple() : bounding_box_complex(rpn_);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
||||
bool
|
||||
|
|
@ -691,6 +792,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
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -983,6 +1094,55 @@ 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) {
|
||||
if (index < 0 || index >= model::cells.size()) {
|
||||
set_errmsg("Index in cells array is out of bounds.");
|
||||
return OPENMC_E_OUT_OF_BOUNDS;
|
||||
}
|
||||
|
||||
*name = model::cells[index]->name().data();
|
||||
|
||||
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()) {
|
||||
set_errmsg("Index in cells array is out of bounds.");
|
||||
return OPENMC_E_OUT_OF_BOUNDS;
|
||||
}
|
||||
|
||||
model::cells[index]->set_name(name);
|
||||
|
||||
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)
|
||||
|
|
|
|||
|
|
@ -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,21 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance)
|
|||
return 0;
|
||||
}
|
||||
|
||||
extern "C" int openmc_global_bounding_box(double* llc, double* urc) {
|
||||
auto bbox = model::universes.at(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
|
||||
|
|
|
|||
|
|
@ -1381,6 +1381,30 @@ 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()) {
|
||||
set_errmsg("Index in materials array is out of bounds.");
|
||||
return OPENMC_E_OUT_OF_BOUNDS;
|
||||
}
|
||||
|
||||
*name = model::materials[index]->name().data();
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
extern "C" int
|
||||
openmc_material_set_name(int32_t index, const char* name) {
|
||||
if (index < 0 || index >= model::materials.size()) {
|
||||
set_errmsg("Index in materials array is out of bounds.");
|
||||
return OPENMC_E_OUT_OF_BOUNDS;
|
||||
}
|
||||
|
||||
model::materials[index]->set_name(name);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
extern "C" int
|
||||
openmc_material_set_volume(int32_t index, double volume)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -271,17 +271,8 @@ Direction DAGSurface::reflect(Position r, Direction u) const
|
|||
return simulation::last_dir;
|
||||
}
|
||||
|
||||
BoundingBox DAGSurface::bounding_box() 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
|
||||
//==============================================================================
|
||||
// PeriodicSurface implementation
|
||||
|
|
@ -366,9 +357,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 +423,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 +473,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_};
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -538,12 +541,6 @@ bool SurfacePlane::periodic_translate(const PeriodicSurface* other, Position& r,
|
|||
return false;
|
||||
}
|
||||
|
||||
BoundingBox
|
||||
SurfacePlane::bounding_box() const
|
||||
{
|
||||
return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Generic functions for x-, y-, and z-, cylinders
|
||||
//==============================================================================
|
||||
|
|
@ -645,7 +642,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 +649,13 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
|
|||
write_dataset(group_id, "coefficients", coeffs);
|
||||
}
|
||||
|
||||
BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const {
|
||||
if (!pos_side) {
|
||||
return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_};
|
||||
} else {
|
||||
return {};
|
||||
}
|
||||
}
|
||||
//==============================================================================
|
||||
// SurfaceYCylinder implementation
|
||||
//==============================================================================
|
||||
|
|
@ -686,6 +689,14 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
|
|||
write_dataset(group_id, "coefficients", coeffs);
|
||||
}
|
||||
|
||||
BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const {
|
||||
if (!pos_side) {
|
||||
return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_, z0_ + radius_};
|
||||
} else {
|
||||
return {};
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// SurfaceZCylinder implementation
|
||||
//==============================================================================
|
||||
|
|
@ -719,6 +730,15 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
|
|||
write_dataset(group_id, "coefficients", coeffs);
|
||||
}
|
||||
|
||||
BoundingBox SurfaceZCylinder::bounding_box(bool pos_side) const {
|
||||
if (!pos_side) {
|
||||
return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY, INFTY};
|
||||
} else {
|
||||
return {};
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
//==============================================================================
|
||||
// SurfaceSphere implementation
|
||||
//==============================================================================
|
||||
|
|
@ -787,6 +807,16 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
|
|||
write_dataset(group_id, "coefficients", coeffs);
|
||||
}
|
||||
|
||||
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 {};
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Generic functions for x-, y-, and z-, cones
|
||||
//==============================================================================
|
||||
|
|
@ -1161,7 +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();
|
||||
auto bb = surf->bounding_box(true) & surf->bounding_box(false);
|
||||
if (bb.xmin > -INFTY && bb.xmin < xmin) {
|
||||
xmin = bb.xmin;
|
||||
i_xmin = i_surf;
|
||||
|
|
|
|||
|
|
@ -1,5 +1,8 @@
|
|||
from tests.testing_harness import TestHarness
|
||||
|
||||
import sys
|
||||
|
||||
import openmc.capi
|
||||
|
||||
def test_complex_cell():
|
||||
harness = TestHarness('statepoint.10.h5')
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -73,7 +73,9 @@ 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"
|
||||
cell.name = "Not fuel"
|
||||
assert cell.name == "Not fuel"
|
||||
|
||||
def test_cell_temperature(capi_init):
|
||||
cell = openmc.capi.cells[1]
|
||||
|
|
@ -122,7 +124,9 @@ 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]
|
||||
|
|
@ -469,3 +473,13 @@ def test_position(capi_init):
|
|||
pos[2] = 3.3
|
||||
|
||||
assert tuple(pos) == (1.3, 2.3, 3.3)
|
||||
|
||||
|
||||
def test_global_bounding_box(capi_init):
|
||||
expected_llc = (-0.63, -0.63, -np.inf)
|
||||
expected_urc = (0.63, 0.63, np.inf)
|
||||
|
||||
llc, urc = openmc.capi.global_bounding_box()
|
||||
|
||||
assert tuple(llc) == expected_llc
|
||||
assert tuple(urc) == expected_urc
|
||||
|
|
|
|||
98
tests/unit_tests/test_complex_cell_capi.py
Normal file
98
tests/unit_tests/test_complex_cell_capi.py
Normal file
|
|
@ -0,0 +1,98 @@
|
|||
import numpy as np
|
||||
import openmc.capi
|
||||
import pytest
|
||||
|
||||
@pytest.fixture(autouse=True)
|
||||
def 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)
|
||||
|
||||
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, c5])
|
||||
|
||||
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()
|
||||
|
||||
yield
|
||||
|
||||
openmc.capi.finalize()
|
||||
|
||||
|
||||
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
|
||||
assert tuple(cell_box[0]) == expected_box[0]
|
||||
assert tuple(cell_box[1]) == expected_box[1]
|
||||
Loading…
Add table
Add a link
Reference in a new issue