Move cell%type to C++

This commit is contained in:
Sterling Harper 2018-05-11 17:20:25 -04:00
parent e19fc61f1d
commit 110ff69cac
12 changed files with 194 additions and 109 deletions

View file

@ -435,6 +435,7 @@ set(LIBOPENMC_CXX_SRC
src/cell.cpp
src/initialize.cpp
src/finalize.cpp
src/geometry_aux.cpp
src/hdf5_interface.cpp
src/lattice.cpp
src/message_passing.cpp

View file

@ -4,7 +4,6 @@
#include <limits>
#include <sstream>
#include <string>
#include <unordered_set>
#include "constants.h"
#include "error.h"
@ -226,6 +225,15 @@ Cell::Cell(pugi::xml_node cell_node)
fill = C_NONE;
}
if (check_for_node(cell_node, "material")) {
//TODO: read material ids.
material.push_back(C_NONE+1);
material.shrink_to_fit();
} else {
material.push_back(C_NONE);
material.shrink_to_fit();
}
std::string region_spec {""};
if (check_for_node(cell_node, "region")) {
region_spec = get_node_value(cell_node, "region");
@ -443,74 +451,6 @@ read_cells(pugi::xml_node *node)
}
}
//==============================================================================
extern "C" int32_t
find_root_universe()
{
// Find all the universes listed as a cell fill.
std::unordered_set<int32_t> fill_univ_ids;
for (Cell *c : cells_c) {
fill_univ_ids.insert(c->fill);
}
// Find all the universes contained in a lattice.
for (Lattice *lat : lattices_c) {
for (auto it = lat->begin(); it != lat->end(); ++it) {
fill_univ_ids.insert(*it);
}
if (lat->outer != NO_OUTER_UNIVERSE) {
fill_univ_ids.insert(lat->outer);
}
}
// Figure out which universe is not in the set. This is the root universe.
bool root_found {false};
int32_t root_univ;
for (int32_t i = 0; i < universes_c.size(); i++) {
auto search = fill_univ_ids.find(universes_c[i]->id);
if (search == fill_univ_ids.end()) {
if (root_found) {
fatal_error("Two or more universes are not used as fill universes, so "
"it is not possible to distinguish which one is the root "
"universe.");
} else {
root_found = true;
root_univ = i;
}
}
}
if (!root_found) fatal_error("Could not find a root universe. Make sure "
"there are no circular dependencies in the geometry.");
return root_univ;
}
//==============================================================================
extern "C" void
adjust_indices_c()
{
// Change cell.universe values from IDs to indices.
for (Cell *c : cells_c) {
auto it = universe_dict.find(c->universe);
if (it != universe_dict.end()) {
//TODO: Remove this off-by-one indexing.
c->universe = it->second + 1;
} else {
std::stringstream err_msg;
err_msg << "Could not find universe " << c->universe
<< " specified on cell " << c->id;
fatal_error(err_msg);
}
}
// Change all lattice universe values from IDs to indices.
for (Lattice *l : lattices_c) {
l->adjust_indices();
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
@ -522,6 +462,10 @@ extern "C" {
void cell_set_id(Cell *c, int32_t id) {c->id = id;}
int cell_type(Cell *c) {return c->type;}
void cell_set_type(Cell *c, int type) {c->type = type;}
int32_t cell_universe(Cell *c) {return c->universe;}
void cell_set_universe(Cell *c, int32_t universe) {c->universe = universe;}

View file

@ -12,6 +12,14 @@
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
extern "C" int FILL_MATERIAL;
extern "C" int FILL_UNIVERSE;
extern "C" int FILL_LATTICE;
//==============================================================================
// Global variables
//==============================================================================
@ -48,9 +56,14 @@ class Cell
public:
int32_t id; //!< Unique ID
std::string name{""}; //!< User-defined name
int type; //!< Material, universe, or lattice
int32_t universe; //!< Universe # this cell is in
int32_t fill; //!< Universe # filling this cell
//! Material within this cell. May be multiple materials for distribcell.
//! C_NONE signifies a universe.
std::vector<int32_t> material;
//! Definition of spatial region as Boolean expression of half-spaces
std::vector<std::int32_t> region;
//! Reverse Polish notation for region expression

View file

@ -113,6 +113,9 @@ module constants
FILL_MATERIAL = 1, & ! Cell with a specified material
FILL_UNIVERSE = 2, & ! Cell filled by a separate universe
FILL_LATTICE = 3 ! Cell filled with a lattice
integer(C_INT), bind(C, name='FILL_MATERIAL') :: FILL_MATERIAL_C = FILL_MATERIAL
integer(C_INT), bind(C, name='FILL_UNIVERSE') :: FILL_UNIVERSE_C = FILL_UNIVERSE
integer(C_INT), bind(C, name='FILL_LATTICE') :: FILL_LATTICE_C = FILL_LATTICE
! Void material
integer, parameter :: MATERIAL_VOID = -1

View file

@ -1,6 +1,8 @@
#ifndef CONSTANTS_H
#define CONSTANTS_H
#include <limits>
namespace openmc{

View file

@ -150,7 +150,7 @@ contains
if (found) then
associate(c => cells(i_cell))
CELL_TYPE: if (c % type == FILL_MATERIAL) then
CELL_TYPE: if (c % type() == FILL_MATERIAL) then
! ======================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
@ -166,10 +166,10 @@ contains
distribcell_index = c % distribcell_index
offset = 0
do k = 1, p % n_coord
if (cells(p % coord(k) % cell) % type == FILL_UNIVERSE) then
if (cells(p % coord(k) % cell) % type() == FILL_UNIVERSE) then
offset = offset + cells(p % coord(k) % cell) % &
offset(distribcell_index)
elseif (cells(p % coord(k) % cell) % type == FILL_LATTICE) then
elseif (cells(p % coord(k) % cell) % type() == FILL_LATTICE) then
if (lattices(p % coord(k + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(k + 1) % lattice_x, &
@ -204,7 +204,7 @@ contains
p % sqrtkT = c % sqrtkT(1)
end if
elseif (c % type == FILL_UNIVERSE) then CELL_TYPE
elseif (c % type() == FILL_UNIVERSE) then CELL_TYPE
! ======================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
@ -232,7 +232,7 @@ contains
call find_cell(p, found)
j = p % n_coord
elseif (c % type == FILL_LATTICE) then CELL_TYPE
elseif (c % type() == FILL_LATTICE) then CELL_TYPE
! ======================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
@ -561,11 +561,11 @@ contains
! ====================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
if (c % type == FILL_MATERIAL) then
if (c % type() == FILL_MATERIAL) then
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
elseif (c % type == FILL_UNIVERSE) then
elseif (c % type() == FILL_UNIVERSE) then
! Set offset for the cell on this level
c % offset(map) = offset
@ -579,7 +579,7 @@ contains
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == FILL_LATTICE) then
elseif (c % type() == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
@ -671,11 +671,11 @@ contains
! ====================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
if (c % type == FILL_MATERIAL) then
if (c % type() == FILL_MATERIAL) then
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
elseif (c % type == FILL_UNIVERSE) then
elseif (c % type() == FILL_UNIVERSE) then
next_univ => universes(c % fill)
@ -690,7 +690,7 @@ contains
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == FILL_LATTICE) then
elseif (c % type() == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
@ -771,13 +771,13 @@ contains
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
if (c % type == FILL_UNIVERSE) then
if (c % type() == FILL_UNIVERSE) then
call count_instance(universes(c % fill))
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == FILL_LATTICE) then
elseif (c % type() == FILL_LATTICE) then
! Set current lattice
associate (lat => lattices(c % fill) % obj)
@ -840,14 +840,14 @@ contains
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
if (c % type == FILL_UNIVERSE) then
if (c % type() == FILL_UNIVERSE) then
next_univ => universes(c % fill)
levels_below = max(levels_below, maximum_levels(next_univ))
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == FILL_LATTICE) then
elseif (c % type() == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj

104
src/geometry_aux.cpp Normal file
View file

@ -0,0 +1,104 @@
#include <sstream>
#include <unordered_set>
#include "cell.h"
#include "constants.h"
#include "error.h"
#include "lattice.h"
#include <iostream> //TODO: remove this
namespace openmc {
extern "C" int32_t
find_root_universe()
{
// Find all the universes listed as a cell fill.
std::unordered_set<int32_t> fill_univ_ids;
for (Cell *c : cells_c) {
fill_univ_ids.insert(c->fill);
}
// Find all the universes contained in a lattice.
for (Lattice *lat : lattices_c) {
for (auto it = lat->begin(); it != lat->end(); ++it) {
fill_univ_ids.insert(*it);
}
if (lat->outer != NO_OUTER_UNIVERSE) {
fill_univ_ids.insert(lat->outer);
}
}
// Figure out which universe is not in the set. This is the root universe.
bool root_found {false};
int32_t root_univ;
for (int32_t i = 0; i < universes_c.size(); i++) {
auto search = fill_univ_ids.find(universes_c[i]->id);
if (search == fill_univ_ids.end()) {
if (root_found) {
fatal_error("Two or more universes are not used as fill universes, so "
"it is not possible to distinguish which one is the root "
"universe.");
} else {
root_found = true;
root_univ = i;
}
}
}
if (!root_found) fatal_error("Could not find a root universe. Make sure "
"there are no circular dependencies in the geometry.");
return root_univ;
}
//==============================================================================
extern "C" void
adjust_indices_c()
{
// Adjust material/fill idices.
for (Cell *c : cells_c) {
if (c->material[0] == C_NONE) {
int32_t id = c->fill;
auto search_univ = universe_dict.find(id);
auto search_lat = lattice_dict.find(id);
if (search_univ != universe_dict.end()) {
c->type = FILL_UNIVERSE;
c->fill = search_univ->second + 1; //TODO: off-by-one
} else if (search_lat != lattice_dict.end()) {
c->type = FILL_LATTICE;
c->fill = search_lat->second + 1; //TODO: off-by-one
} else {
std::stringstream err_msg;
err_msg << "Specified fill " << id << " on cell " << c->id
<< " is neither a universe nor a lattice.";
fatal_error(err_msg);
}
} else {
//TODO: materials
c->type = FILL_MATERIAL;
}
}
// Change cell.universe values from IDs to indices.
for (Cell *c : cells_c) {
auto search = universe_dict.find(c->universe);
if (search != universe_dict.end()) {
//TODO: Remove this off-by-one indexing.
c->universe = search->second + 1;
} else {
std::stringstream err_msg;
err_msg << "Could not find universe " << c->universe
<< " specified on cell " << c->id;
fatal_error(err_msg);
}
}
// Change all lattice universe values from IDs to indices.
for (Lattice *l : lattices_c) {
l->adjust_indices();
}
}
} // namespace openmc

View file

@ -37,6 +37,20 @@ module geometry_header
integer(C_INT32_T), intent(in), value :: id
end subroutine cell_set_id_c
function cell_type_c(cell_ptr) bind(C, name='cell_type') result(type)
import C_PTR, C_INT
implicit none
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT) :: type
end function cell_type_c
subroutine cell_set_type_c(cell_ptr, type) bind(C, name='cell_set_type')
import C_PTR, C_INT
implicit none
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT), intent(in), value :: type
end subroutine cell_set_type_c
function cell_universe_c(cell_ptr) bind(C, name='cell_universe') &
result(universe)
import C_PTR, C_INT32_T
@ -230,8 +244,6 @@ module geometry_header
type Cell
type(C_PTR) :: ptr
integer :: type ! Type of cell (normal, universe,
! lattice)
integer :: fill ! universe # filling this cell
integer :: instances ! number of instances of this cell in
! the geom
@ -257,6 +269,8 @@ module geometry_header
procedure :: id => cell_id
procedure :: set_id => cell_set_id
procedure :: type => cell_type
procedure :: set_type => cell_set_type
procedure :: universe => cell_universe
procedure :: set_universe => cell_set_universe
procedure :: simple => cell_simple
@ -355,6 +369,18 @@ contains
call cell_set_id_c(this % ptr, id)
end subroutine cell_set_id
function cell_type(this) result(type)
class(Cell), intent(in) :: this
integer(C_INT) :: type
type = cell_type_c(this % ptr)
end function cell_type
subroutine cell_set_type(this, type)
class(Cell), intent(in) :: this
integer(C_INT), intent(in) :: type
call cell_set_type_c(this % ptr, type)
end subroutine cell_set_type
function cell_universe(this) result(universe)
class(Cell), intent(in) :: this
integer(C_INT32_T) :: universe
@ -544,7 +570,7 @@ contains
err = 0
if (index >= 1 .and. index <= size(cells)) then
associate (c => cells(index))
type = c % type
type = c % type()
select case (type)
case (FILL_MATERIAL)
n = size(c % material)
@ -595,7 +621,7 @@ contains
if (allocated(c % material)) deallocate(c % material)
allocate(c % material(n))
c % type = FILL_MATERIAL
call c % set_type(FILL_MATERIAL)
do i = 1, n
j = indices(i)
if ((j >= 1 .and. j <= n_materials) .or. j == MATERIAL_VOID) then
@ -607,9 +633,9 @@ contains
end if
end do
case (FILL_UNIVERSE)
c % type = FILL_UNIVERSE
call c % set_type(FILL_UNIVERSE)
case (FILL_LATTICE)
c % type = FILL_LATTICE
call c % set_type(FILL_LATTICE)
end select
end associate
else

View file

@ -3794,24 +3794,16 @@ contains
if (c % material(1) == NONE) then
id = c % fill
if (universe_dict % has(id)) then
c % type = FILL_UNIVERSE
c % fill = universe_dict % get(id)
elseif (lattice_dict % has(id)) then
lid = lattice_dict % get(id)
c % type = FILL_LATTICE
c % fill = lid
else
call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "&
// trim(to_str(c % id())) // " is neither a universe nor a &
&lattice.")
end if
else
do j = 1, size(c % material)
id = c % material(j)
if (id == MATERIAL_VOID) then
c % type = FILL_MATERIAL
else if (material_dict % has(id)) then
c % type = FILL_MATERIAL
c % material(j) = material_dict % get(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
@ -3989,7 +3981,7 @@ contains
! Allocate offset table for fill cells
do i = 1, n_cells
if (cells(i) % type /= FILL_MATERIAL) then
if (cells(i) % type() /= FILL_MATERIAL) then
allocate(cells(i) % offset(n_maps))
end if
end do

View file

@ -85,7 +85,7 @@ contains
if (pl % color_by == PLOT_COLOR_MATS) then
! Assign color based on material
associate (c => cells(p % coord(j) % cell))
if (c % type == FILL_UNIVERSE) then
if (c % type() == FILL_UNIVERSE) then
! If we stopped on a middle universe level, treat as if not found
rgb = pl % not_found % rgb
id = -1

View file

@ -193,7 +193,7 @@ contains
call c % to_hdf5(cell_group)
! Write information on what fills this cell
select case (c%type)
select case (c%type())
case (FILL_MATERIAL)
call write_dataset(cell_group, "fill_type", "material")

View file

@ -58,10 +58,10 @@ contains
distribcell_index = cells(this % cell) % distribcell_index
offset = 0
do i = 1, p % n_coord
if (cells(p % coord(i) % cell) % type == FILL_UNIVERSE) then
if (cells(p % coord(i) % cell) % type() == FILL_UNIVERSE) then
offset = offset + cells(p % coord(i) % cell) % &
offset(distribcell_index)
elseif (cells(p % coord(i) % cell) % type == FILL_LATTICE) then
elseif (cells(p % coord(i) % cell) % type() == FILL_LATTICE) then
if (lattices(p % coord(i + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(i + 1) % lattice_x, &
@ -197,20 +197,20 @@ contains
c => cells(univ % cells(j))
! Skip normal cells which do not have offsets
if (c % type == FILL_MATERIAL) cycle
if (c % type() == FILL_MATERIAL) cycle
! Break loop once we've found the next cell with an offset
exit
end do
! Ensure we didn't just end the loop by iteration
if (c % type /= FILL_MATERIAL) then
if (c % type() /= FILL_MATERIAL) then
! There are more cells in this universe that it could be in
later_cell = .true.
! Two cases, lattice or fill cell
if (c % type == FILL_UNIVERSE) then
if (c % type() == FILL_UNIVERSE) then
temp_offset = c % offset(map)
! Get the offset of the first lattice location
@ -227,7 +227,7 @@ contains
end if
end if
if (n == 1 .and. c % type /= FILL_MATERIAL) then
if (n == 1 .and. c % type() /= FILL_MATERIAL) then
this_cell = .true.
end if
@ -245,7 +245,7 @@ contains
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
if (c % type == FILL_UNIVERSE) then
if (c % type() == FILL_UNIVERSE) then
! Enter this cell to update the current offset
offset = c % offset(map) + offset
@ -256,7 +256,7 @@ contains
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == FILL_LATTICE) then
elseif (c % type() == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj