#include "openmc/cell.h"

#include <algorithm>

#include <cassert>

#include <cctype>

#include <cmath>

#include <iterator>

#include <set>

#include <sstream>

#include <string>

#include <fmt/core.h>

#include "openmc/capi.h"

#include "openmc/constants.h"

#include "openmc/dagmc.h"

#include "openmc/error.h"

#include "openmc/geometry.h"

#include "openmc/hdf5_interface.h"

#include "openmc/lattice.h"

#include "openmc/material.h"

#include "openmc/nuclide.h"

#include "openmc/settings.h"

#include "openmc/xml_interface.h"

namespace openmc {

//==============================================================================

// Global variables

//==============================================================================

namespace model {

std::unordered_map<int32_t, int32_t> cell_map;

vector<unique_ptr<Cell>> cells;

} // namespace model

//==============================================================================

// Cell implementation

//==============================================================================

int32_t Cell::n_instances() const

{

return model::universes[universe_]->n_instances_;

}

void Cell::set_rotation(const vector<double>& rot)

{

if (fill_ == C_NONE) {

fatal_error(fmt::format("Cannot apply a rotation to cell {}"

" because it is not filled with another universe",

id_));

}

if (rot.size() != 3 && rot.size() != 9) {

fatal_error(fmt::format("Non-3D rotation vector applied to cell {}", id_));

}

// Compute and store the inverse rotation matrix for the angles given.

rotation_.clear();

rotation_.reserve(rot.size() == 9 ? 9 : 12);

if (rot.size() == 3) {

double phi = -rot[0] * PI / 180.0;

double theta = -rot[1] * PI / 180.0;

double psi = -rot[2] * PI / 180.0;

rotation_.push_back(std::cos(theta) * std::cos(psi));

rotation_.push_back(-std::cos(phi) * std::sin(psi) +

std::sin(phi) * std::sin(theta) * std::cos(psi));

rotation_.push_back(std::sin(phi) * std::sin(psi) +

std::cos(phi) * std::sin(theta) * std::cos(psi));

rotation_.push_back(std::cos(theta) * std::sin(psi));

rotation_.push_back(std::cos(phi) * std::cos(psi) +

std::sin(phi) * std::sin(theta) * std::sin(psi));

rotation_.push_back(-std::sin(phi) * std::cos(psi) +

std::cos(phi) * std::sin(theta) * std::sin(psi));

rotation_.push_back(-std::sin(theta));

rotation_.push_back(std::sin(phi) * std::cos(theta));

rotation_.push_back(std::cos(phi) * std::cos(theta));

// When user specifies angles, write them at end of vector

rotation_.push_back(rot[0]);

rotation_.push_back(rot[1]);

rotation_.push_back(rot[2]);

} else {

std::copy(rot.begin(), rot.end(), std::back_inserter(rotation_));

}

}

double Cell::temperature(int32_t instance) const

{

if (sqrtkT_.size() < 1) {

throw std::runtime_error {"Cell temperature has not yet been set."};

}

if (instance >= 0) {

double sqrtkT = sqrtkT_.size() == 1 ? sqrtkT_.at(0) : sqrtkT_.at(instance);

return sqrtkT * sqrtkT / K_BOLTZMANN;

} else {

return sqrtkT_[0] * sqrtkT_[0] / K_BOLTZMANN;

}

}

double Cell::density_mult(int32_t instance) const

{

if (instance >= 0) {

return density_mult_.size() == 1 ? density_mult_.at(0)

: density_mult_.at(instance);

} else {

return density_mult_[0];

}

}

double Cell::density(int32_t instance) const

{

const int32_t mat_index = material(instance);

if (mat_index == MATERIAL_VOID)

return 0.0;

return density_mult(instance) * model::materials[mat_index]->density_gpcc();

}

void Cell::set_temperature(double T, int32_t instance, bool set_contained)

{

if (settings::temperature_method == TemperatureMethod::INTERPOLATION) {

if (T < (data::temperature_min - settings::temperature_tolerance)) {

throw std::runtime_error {

fmt::format("Temperature of {} K is below minimum temperature at "

"which data is available of {} K.",

T, data::temperature_min)};

} else if (T > (data::temperature_max + settings::temperature_tolerance)) {

throw std::runtime_error {

fmt::format("Temperature of {} K is above maximum temperature at "

"which data is available of {} K.",

T, data::temperature_max)};

}

}

if (type_ == Fill::MATERIAL) {

if (instance >= 0) {

// If temperature vector is not big enough, resize it first

if (sqrtkT_.size() != n_instances())

sqrtkT_.resize(n_instances(), sqrtkT_[0]);

// Set temperature for the corresponding instance

sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T);

} else {

// Set temperature for all instances

for (auto& T_ : sqrtkT_) {

T_ = std::sqrt(K_BOLTZMANN * T);

}

}

} else {

if (!set_contained) {

throw std::runtime_error {

fmt::format("Attempted to set the temperature of cell {} "

"which is not filled by a material.",

id_)};

}

auto contained_cells = this->get_contained_cells(instance);

for (const auto& entry : contained_cells) {

auto& cell = model::cells[entry.first];

assert(cell->type_ == Fill::MATERIAL);

auto& instances = entry.second;

for (auto instance : instances) {

cell->set_temperature(T, instance);

}

}

}

}

void Cell::set_density(double density, int32_t instance, bool set_contained)

{

if (type_ != Fill::MATERIAL && !set_contained) {

fatal_error(

fmt::format("Attempted to set the density multiplier of cell {} "

"which is not filled by a material.",

id_));

}

if (type_ == Fill::MATERIAL) {

const int32_t mat_index = material(instance);

if (mat_index == MATERIAL_VOID)

return;

if (instance >= 0) {

// If density multiplier vector is not big enough, resize it first

if (density_mult_.size() != n_instances())

density_mult_.resize(n_instances(), density_mult_[0]);

// Set density multiplier for the corresponding instance

density_mult_.at(instance) =

density / model::materials[mat_index]->density_gpcc();

} else {

// Set density multiplier for all instances

for (auto& x : density_mult_) {

x = density / model::materials[mat_index]->density_gpcc();

}

}

} else {

auto contained_cells = this->get_contained_cells(instance);

for (const auto& entry : contained_cells) {

auto& cell = model::cells[entry.first];

assert(cell->type_ == Fill::MATERIAL);

auto& instances = entry.second;

for (auto instance : instances) {

cell->set_density(density, instance);

}

}

}

}

void Cell::export_properties_hdf5(hid_t group) const

{

// Create a group for this cell.

auto cell_group = create_group(group, fmt::format("cell {}", id_));

// Write temperature in [K] for one or more cell instances

vector<double> temps;

for (auto sqrtkT_val : sqrtkT_)

temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN);

write_dataset(cell_group, "temperature", temps);

// Write density for one or more cell instances

if (type_ == Fill::MATERIAL && material_.size() > 0) {

vector<double> density;

for (int32_t i = 0; i < density_mult_.size(); ++i)

density.push_back(this->density(i));

write_dataset(cell_group, "density", density);

}

close_group(cell_group);

}

void Cell::import_properties_hdf5(hid_t group)

{

auto cell_group = open_group(group, fmt::format("cell {}", id_));

// Read temperatures from file

vector<double> temps;

read_dataset(cell_group, "temperature", temps);

// Ensure number of temperatures makes sense

auto n_temps = temps.size();

if (n_temps > 1 && n_temps != n_instances()) {

fatal_error(fmt::format(

"Number of temperatures for cell {} doesn't match number of instances",

id_));

}

// Modify temperatures for the cell

sqrtkT_.clear();

sqrtkT_.resize(temps.size());

for (int64_t i = 0; i < temps.size(); ++i) {

this->set_temperature(temps[i], i);

}

// Read densities

if (object_exists(cell_group, "density")) {

vector<double> density;

read_dataset(cell_group, "density", density);

// Ensure number of densities makes sense

auto n_density = density.size();

if (n_density > 1 && n_density != n_instances()) {

fatal_error(fmt::format("Number of densities for cell {} "

"doesn't match number of instances",

id_));

}

// Set densities.

for (int32_t i = 0; i < n_density; ++i) {

this->set_density(density[i], i);

}

}

close_group(cell_group);

}

void Cell::to_hdf5(hid_t cell_group) const

{

// Create a group for this cell.

auto group = create_group(cell_group, fmt::format("cell {}", id_));

if (!name_.empty()) {

write_string(group, "name", name_, false);

}

write_dataset(group, "universe", model::universes[universe_]->id_);

to_hdf5_inner(group);

// Write fill information.

if (type_ == Fill::MATERIAL) {

write_dataset(group, "fill_type", "material");

std::vector<int32_t> mat_ids;

for (auto i_mat : material_) {

if (i_mat != MATERIAL_VOID) {

mat_ids.push_back(model::materials[i_mat]->id_);

} else {

mat_ids.push_back(MATERIAL_VOID);

}

}

if (mat_ids.size() == 1) {

write_dataset(group, "material", mat_ids[0]);

} else {

write_dataset(group, "material", mat_ids);

}

std::vector<double> temps;

for (auto sqrtkT_val : sqrtkT_)

temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN);

write_dataset(group, "temperature", temps);

write_dataset(group, "density_mult", density_mult_);

} else if (type_ == Fill::UNIVERSE) {

write_dataset(group, "fill_type", "universe");

write_dataset(group, "fill", model::universes[fill_]->id_);

if (translation_ != Position(0, 0, 0)) {

write_dataset(group, "translation", translation_);

}

if (!rotation_.empty()) {

if (rotation_.size() == 12) {

std::array<double, 3> rot {rotation_[9], rotation_[10], rotation_[11]};

write_dataset(group, "rotation", rot);

} else {

write_dataset(group, "rotation", rotation_);

}

}

} else if (type_ == Fill::LATTICE) {

write_dataset(group, "fill_type", "lattice");

write_dataset(group, "lattice", model::lattices[fill_]->id_);

}

close_group(group);

}

//==============================================================================

// XML parsing helpers for <cell> nodes

//==============================================================================

vector<int32_t> parse_cell_material_xml(pugi::xml_node node, int32_t cell_id)

{

vector<std::string> mats {

get_node_array<std::string>(node, "material", true)};

if (mats.empty()) {

fatal_error(fmt::format(

"An empty material element was specified for cell {}", cell_id));

}

vector<int32_t> material;

material.reserve(mats.size());

for (const auto& mat : mats) {

if (mat == "void") {

material.push_back(MATERIAL_VOID);

} else {

material.push_back(std::stoi(mat));

}

}

return material;

}

vector<double> parse_cell_temperature_xml(pugi::xml_node node, int32_t cell_id)

{

auto temperatures = get_node_array<double>(node, "temperature");

if (temperatures.empty()) {

fatal_error(fmt::format(

"An empty temperature element was specified for cell {}", cell_id));

}

for (auto T : temperatures) {

if (T < 0) {

fatal_error(fmt::format(

"Cell {} was specified with a negative temperature", cell_id));

}

}

return temperatures;

}

vector<double> parse_cell_density_xml(pugi::xml_node node, int32_t cell_id)

{

auto densities = get_node_array<double>(node, "density");

if (densities.empty()) {

fatal_error(fmt::format(

"An empty density element was specified for cell {}", cell_id));

}

for (auto rho : densities) {

if (rho <= 0) {

fatal_error(fmt::format(

"Cell {} was specified with a density less than or equal to zero",

cell_id));

}

}

return densities;

}

//==============================================================================

// CSGCell implementation

//==============================================================================

CSGCell::CSGCell(pugi::xml_node cell_node)

{

if (check_for_node(cell_node, "id")) {

id_ = std::stoi(get_node_value(cell_node, "id"));

} else {

fatal_error("Must specify id of cell in geometry XML file.");

}

if (check_for_node(cell_node, "name")) {

name_ = get_node_value(cell_node, "name");

}

if (check_for_node(cell_node, "universe")) {

universe_ = std::stoi(get_node_value(cell_node, "universe"));

} else {

universe_ = 0;

}

// Make sure that either material or fill was specified, but not both.

bool fill_present = check_for_node(cell_node, "fill");

bool material_present = check_for_node(cell_node, "material");

if (!(fill_present || material_present)) {

fatal_error(

fmt::format("Neither material nor fill was specified for cell {}", id_));

}

if (fill_present && material_present) {

fatal_error(fmt::format("Cell {} has both a material and a fill specified; "

"only one can be specified per cell",

id_));

}

if (fill_present) {

fill_ = std::stoi(get_node_value(cell_node, "fill"));

if (fill_ == universe_) {

fatal_error(fmt::format("Cell {} is filled with the same universe that "

"it is contained in.",

id_));

}

} else {

fill_ = C_NONE;

}

// Read the material element. There can be zero materials (filled with a

// universe), more than one material (distribmats), and some materials may

// be "void".

if (material_present) {

material_ = parse_cell_material_xml(cell_node, id_);

}

// Read the temperature element which may be distributed like materials.

if (check_for_node(cell_node, "temperature")) {

sqrtkT_ = parse_cell_temperature_xml(cell_node, id_);

sqrtkT_.shrink_to_fit();

// Make sure this is a material-filled cell.

if (material_.size() == 0) {

fatal_error(fmt::format(

"Cell {} was specified with a temperature but no material. Temperature"

"specification is only valid for cells filled with a material.",

id_));

}

// Convert to sqrt(k*T).

for (auto& T : sqrtkT_) {

T = std::sqrt(K_BOLTZMANN * T);

}

}

// Read the density element which can be distributed similar to temperature.

// These get assigned to the density multiplier, requiring a division by

// the material density.

// Note: calculating the actual density multiplier is deferred until materials

// are finalized. density_mult_ contains the true density in the meantime.

if (check_for_node(cell_node, "density")) {

density_mult_ = parse_cell_density_xml(cell_node, id_);

density_mult_.shrink_to_fit();

// Make sure this is a material-filled cell.

if (material_.size() == 0) {

fatal_error(fmt::format(

"Cell {} was specified with a density but no material. Density"

"specification is only valid for cells filled with a material.",

id_));

}

// Make sure this is a non-void material.

for (auto mat_id : material_) {

if (mat_id == MATERIAL_VOID) {

fatal_error(fmt::format(

"Cell {} was specified with a density, but contains a void "

"material. Density specification is only valid for cells "

"filled with a non-void material.",

id_));

}

}

}

// Read the region specification.

std::string region_spec;

if (check_for_node(cell_node, "region")) {

region_spec = get_node_value(cell_node, "region");

}

// Get a tokenized representation of the region specification and apply De

// Morgans law

Region region(region_spec, id_);

region_ = region;

// Read the translation vector.

if (check_for_node(cell_node, "translation")) {

if (fill_ == C_NONE) {

fatal_error(fmt::format("Cannot apply a translation to cell {}"

" because it is not filled with another universe",

id_));

}

auto xyz {get_node_array<double>(cell_node, "translation")};

if (xyz.size() != 3) {

fatal_error(

fmt::format("Non-3D translation vector applied to cell {}", id_));

}

translation_ = xyz;

}

// Read the rotation transform.

if (check_for_node(cell_node, "rotation")) {

auto rot {get_node_array<double>(cell_node, "rotation")};

set_rotation(rot);

}

}

//==============================================================================

void CSGCell::to_hdf5_inner(hid_t group_id) const

{

write_string(group_id, "geom_type", "csg", false);

write_string(group_id, "region", region_.str(), false);

}

//==============================================================================

vector<int32_t>::iterator CSGCell::find_left_parenthesis(

vector<int32_t>::iterator start, const vector<int32_t>& infix)

{

// start search at zero

int parenthesis_level = 0;

auto it = start;

while (it != infix.begin()) {

// look at two tokens at a time

int32_t one = *it;

int32_t two = *(it - 1);

// decrement parenthesis level if there are two adjacent surfaces

if (one < OP_UNION && two < OP_UNION) {

parenthesis_level--;

// increment if there are two adjacent operators

} else if (one >= OP_UNION && two >= OP_UNION) {

parenthesis_level++;

}

// if the level gets to zero, return the position

if (parenthesis_level == 0) {

// move the iterator back one before leaving the loop

// so that all tokens in the parenthesis block are included

it--;

break;

}

// continue loop, one token at a time

it--;

}

return it;

}

//==============================================================================

// Region implementation

//==============================================================================

Region::Region(std::string region_spec, int32_t cell_id)

{

// Check if region_spec is not empty.

if (!region_spec.empty()) {

// Parse all halfspaces and operators except for intersection (whitespace).

for (int i = 0; i < region_spec.size();) {

if (region_spec[i] == '(') {

expression_.push_back(OP_LEFT_PAREN);

i++;

} else if (region_spec[i] == ')') {

expression_.push_back(OP_RIGHT_PAREN);

i++;

} else if (region_spec[i] == '|') {

expression_.push_back(OP_UNION);

i++;

} else if (region_spec[i] == '~') {

expression_.push_back(OP_COMPLEMENT);

i++;

} else if (region_spec[i] == '-' || region_spec[i] == '+' ||

std::isdigit(region_spec[i])) {

// This is the start of a halfspace specification. Iterate j until we

// find the end, then push-back everything between i and j.

int j = i + 1;

while (j < region_spec.size() && std::isdigit(region_spec[j])) {

j++;

}

expression_.push_back(std::stoi(region_spec.substr(i, j - i)));

i = j;

} else if (std::isspace(region_spec[i])) {

i++;

} else {

auto err_msg =

fmt::format("Region specification contains invalid character, \"{}\"",

region_spec[i]);

fatal_error(err_msg);

}

}

// Add in intersection operators where a missing operator is needed.

int i = 0;

while (i < expression_.size() - 1) {

bool left_compat {

(expression_[i] < OP_UNION) || (expression_[i] == OP_RIGHT_PAREN)};

bool right_compat {(expression_[i + 1] < OP_UNION) ||

(expression_[i + 1] == OP_LEFT_PAREN) ||

(expression_[i + 1] == OP_COMPLEMENT)};

if (left_compat && right_compat) {

expression_.insert(expression_.begin() + i + 1, OP_INTERSECTION);

}

i++;

}

// Remove complement operators using DeMorgan's laws

auto it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT);

while (it != expression_.end()) {

// Erase complement. Note that erase invalidates the iterator, so we have

// to use the iterator it returns, which points to the token that

// followed the complement operator.

it = expression_.erase(it);

if (it == expression_.end())

break;

// Define stop given left parenthesis or not

auto stop = it;

if (*it == OP_LEFT_PAREN) {

int depth = 1;

do {

stop++;

if (*stop > OP_COMPLEMENT) {

if (*stop == OP_RIGHT_PAREN) {

depth--;

} else {

depth++;

}

}

} while (depth > 0);

it++;

}

// apply DeMorgan's law to any surfaces/operators between these

// positions in the RPN

apply_demorgan(it, stop);

// update iterator position

it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT);

}

// Convert user IDs to surface indices.

for (auto& r : expression_) {

if (r < OP_UNION) {

const auto& it {model::surface_map.find(abs(r))};

if (it == model::surface_map.end()) {

throw std::runtime_error {

"Invalid surface ID " + std::to_string(abs(r)) +

" specified in region for cell " + std::to_string(cell_id) + "."};

}

r = (r > 0) ? it->second + 1 : -(it->second + 1);

}

}

// Check if this is a simple cell.

simple_ = true;

for (int32_t token : expression_) {

if (token == OP_UNION) {

simple_ = false;

// Ensure intersections have precedence over unions

enforce_precedence();

break;

}

}

// If this cell is simple, remove all the superfluous operator tokens.

if (simple_) {

expression_.erase(std::remove_if(expression_.begin(), expression_.end(),

[](int32_t token) {

return token == OP_INTERSECTION ||

token > OP_COMPLEMENT;

}),

expression_.end());

}

expression_.shrink_to_fit();

} else {

simple_ = true;

}

}

//==============================================================================

void Region::apply_demorgan(

vector<int32_t>::iterator start, vector<int32_t>::iterator stop)

{

do {

if (*start < OP_UNION) {

*start *= -1;

} else if (*start == OP_UNION) {

*start = OP_INTERSECTION;

} else if (*start == OP_INTERSECTION) {

*start = OP_UNION;

}

start++;

} while (start < stop);

}

//==============================================================================

//! Add precedence for infix regions so intersections have higher

//! precedence than unions using parentheses.

//==============================================================================

void Region::add_parentheses(int64_t start)

{

int32_t start_token = expression_[start];

// Add left parenthesis and set new position to be after parenthesis

if (start_token == OP_UNION) {

start += 2;

}

expression_.insert(expression_.begin() + start - 1, OP_LEFT_PAREN);

// Add right parenthesis

// While the start iterator is within the bounds of infix

while (start + 1 < expression_.size()) {

start++;

// If the current token is an operator and is different than the start token

if (expression_[start] >= OP_UNION && expression_[start] != start_token) {

// Skip wrapped regions but save iterator position to check precedence and

// add right parenthesis, right parenthesis position depends on the

// operator, when the operator is a union then do not include the operator

// in the region, when the operator is an intersection then include the

// operator and next surface

if (expression_[start] == OP_LEFT_PAREN) {

int depth = 1;

do {

start++;

if (expression_[start] > OP_COMPLEMENT) {

if (expression_[start] == OP_RIGHT_PAREN) {

depth--;

} else {

depth++;

}

}

} while (depth > 0);

} else {

if (start_token == OP_UNION) {

--start;

}

expression_.insert(expression_.begin() + start, OP_RIGHT_PAREN);

return;

}

}

}

// If we get here a right parenthesis hasn't been placed

expression_.push_back(OP_RIGHT_PAREN);

}

//==============================================================================

//! Add parentheses to enforce operator precedence in region expressions

//!

//! This function ensures that intersection operators have higher precedence

//! than union operators by adding parentheses where needed. For example:

//! "1 2 | 3" becomes "(1 2) | 3"

//! "1 | 2 3" becomes "1 | (2 3)"

//!

//! The algorithm uses stacks to track the current operator type and its

//! position at each parenthesis depth level. When it encounters a different

//! operator at the same depth, it adds parentheses to group the

//! higher-precedence operations.

//==============================================================================

void Region::enforce_precedence()

{

// Stack tracking the operator type at each depth (0 = no operator seen yet)

vector<int32_t> op_stack = {0};

// Stack tracking where the operator sequence started at each depth

vector<std::size_t> pos_stack = {0};

for (int64_t i = 0; i < expression_.size(); ++i) {

int32_t token = expression_[i];

if (token == OP_LEFT_PAREN) {

// Entering a new parenthesis level - push new tracking state

op_stack.push_back(0);

pos_stack.push_back(0);

continue;

} else if (token == OP_RIGHT_PAREN) {

// Exiting a parenthesis level - pop tracking state (keep at least one)

if (op_stack.size() > 1) {

op_stack.pop_back();

pos_stack.pop_back();

}

continue;

}

if (token == OP_UNION || token == OP_INTERSECTION) {

if (op_stack.back() == 0) {

// First operator at this depth - record it and its position

op_stack.back() = token;

pos_stack.back() = i;

} else if (token != op_stack.back()) {

// Encountered a different operator at the same depth - need to add

// parentheses to enforce precedence. Intersection has higher

// precedence, so we parenthesize the intersection terms.

if (op_stack.back() == OP_INTERSECTION) {

add_parentheses(pos_stack.back());

} else {

add_parentheses(i);

}

// Restart the scan since we modified the expression

i = -1; // Will be incremented to 0 by the for loop

op_stack = {0};

pos_stack = {0};

}

}

}

}

//==============================================================================

//! Convert infix region specification to Reverse Polish Notation (RPN)

//!

//! This function uses the shunting-yard algorithm.

//==============================================================================

vector<int32_t> Region::generate_postfix(int32_t cell_id) const

{

vector<int32_t> rpn;

vector<int32_t> stack;

for (int32_t token : expression_) {

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) {

int32_t op = stack.back();

if (op < OP_RIGHT_PAREN && ((token == OP_COMPLEMENT && token < op) ||

(token != OP_COMPLEMENT && token <= op))) {

// While there is an operator, op, on top of the stack, if the token

// is left-associative and its precedence is less than or equal to

// that of op or if the token is right-associative and its precedence

// 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);

stack.pop_back();

} else {

break;

}

}

stack.push_back(token);

} else if (token == OP_LEFT_PAREN) {

// If the token is a left parenthesis, push it onto the stack

stack.push_back(token);

} else {

// If the token is a right parenthesis, move operators from the stack to

// the output queue until reaching the left parenthesis.

for (auto it = stack.rbegin(); *it != OP_LEFT_PAREN; it++) {

// If we run out of operators without finding a left parenthesis, it

// means there are mismatched parentheses.

if (it == stack.rend()) {

fatal_error(fmt::format(

"Mismatched parentheses in region specification for cell {}",

cell_id));

}

rpn.push_back(stack.back());

stack.pop_back();

}

// Pop the left parenthesis.

stack.pop_back();

}

}

while (stack.size() > 0) {

int32_t op = stack.back();

// If the operator is a parenthesis it is mismatched.

if (op >= OP_RIGHT_PAREN) {

fatal_error(fmt::format(

"Mismatched parentheses in region specification for cell {}", cell_id));

}

rpn.push_back(stack.back());

stack.pop_back();

}

return rpn;

}

//==============================================================================

std::string Region::str() const

{

std::stringstream region_spec {};

if (!expression_.empty()) {

for (int32_t token : expression_) {

if (token == OP_LEFT_PAREN) {

region_spec << " (";

} else if (token == OP_RIGHT_PAREN) {

region_spec << " )";

} else if (token == OP_COMPLEMENT) {

region_spec << " ~";

} else if (token == OP_INTERSECTION) {

} else if (token == OP_UNION) {

region_spec << " |";

} else {

// Note the off-by-one indexing

auto surf_id = model::surfaces[abs(token) - 1]->id_;

region_spec << " " << ((token > 0) ? surf_id : -surf_id);

}

}

}

return region_spec.str();

}

//==============================================================================

std::pair<double, int32_t> Region::distance(

Position r, Direction u, int32_t on_surface) const

{

if (simple_) {

return distance_to_nearest_surface(r, u, on_surface, false);

} else {

return distance_complex(r, u, on_surface);

}

}

//==============================================================================

std::pair<double, int32_t> Region::distance_to_nearest_surface(Position r,

Direction u, int32_t on_surface, bool ignore_coincident_surfaces) const

{

double min_dist {INFTY};

int32_t i_surf {std::numeric_limits<int32_t>::max()};

for (int32_t token : expression_) {

// Ignore this token if it corresponds to an operator rather than a region.

if (token >= OP_UNION)

continue;

// Calculate the distance to this surface.

// Note the off-by-one indexing

bool coincident {std::abs(token) == std::abs(on_surface)};

double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};

// Different surface definitions can represent the same geometric surface.

// When the ray is already known to be on a surface, ignore intersections

// with other surfaces at the same location to avoid repeatedly crossing

// between them due to roundoff.

if (ignore_coincident_surfaces && d < FP_COINCIDENT)

continue;

// Check if this distance is the new minimum.

if (d < min_dist) {

if (min_dist - d >= FP_PRECISION * min_dist) {

min_dist = d;

i_surf = -token;

}

}

}

return {min_dist, i_surf};

}

//==============================================================================

std::pair<double, int32_t> Region::distance_complex(

Position r, Direction u, int32_t on_surface) const