Merge pull request #1045 from smharper/cpp_materials

Add Material objects to C++
This commit is contained in:
Paul Romano 2018-08-17 09:02:46 -05:00 committed by GitHub
commit e4232536df
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
23 changed files with 572 additions and 450 deletions

View file

@ -397,6 +397,7 @@ add_library(libopenmc SHARED
src/geometry_aux.cpp
src/hdf5_interface.cpp
src/lattice.cpp
src/material.cpp
src/math_functions.cpp
src/message_passing.cpp
src/mgxs.cpp

View file

@ -106,9 +106,13 @@ class Cell(_FortranObjectWithID):
if fill_type.value == 1:
if n.value > 1:
return [Material(index=i) for i in indices[:n.value]]
#TODO: off-by-one
return [Material(index=i+1 if i >= 0 else i)
for i in indices[:n.value]]
else:
return Material(index=indices[0])
#TODO: off-by-one
index = indices[0] + 1 if indices[0] >= 0 else indices[0]
return Material(index=index)
else:
raise NotImplementedError

View file

@ -37,8 +37,6 @@ module openmc_api
public :: openmc_calculate_volumes
public :: openmc_cell_filter_get_bins
public :: openmc_cell_get_id
public :: openmc_cell_get_fill
public :: openmc_cell_set_fill
public :: openmc_cell_set_id
public :: openmc_cell_set_temperature
public :: openmc_energy_filter_get_bins
@ -214,7 +212,7 @@ contains
if (p % material == MATERIAL_VOID) then
id = 0
else
id = materials(p % material) % id
id = materials(p % material) % id()
end if
end if
instance = p % cell_instance - 1

View file

@ -8,6 +8,8 @@
#include "error.h"
#include "hdf5_interface.h"
#include "lattice.h"
#include "material.h"
#include "openmc.h"
#include "surface.h"
#include "xml_interface.h"
@ -197,52 +199,64 @@ generate_rpn(int32_t cell_id, std::vector<int32_t> infix)
Cell::Cell(pugi::xml_node cell_node)
{
if (check_for_node(cell_node, "id")) {
id = stoi(get_node_value(cell_node, "id"));
id = std::stoi(get_node_value(cell_node, "id"));
} else {
fatal_error("Must specify id of cell in geometry XML file.");
}
//TODO: don't automatically lowercase cell and surface names
if (check_for_node(cell_node, "name")) {
name = get_node_value(cell_node, "name");
}
if (check_for_node(cell_node, "universe")) {
universe = stoi(get_node_value(cell_node, "universe"));
universe = std::stoi(get_node_value(cell_node, "universe"));
} else {
universe = 0;
}
if (check_for_node(cell_node, "fill")) {
fill = stoi(get_node_value(cell_node, "fill"));
} else {
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();
}
// Make sure that either material or fill was specified.
if ((material[0] == C_NONE) && (fill == C_NONE)) {
// 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)) {
std::stringstream err_msg;
err_msg << "Neither material nor fill was specified for cell " << id;
fatal_error(err_msg);
}
// Make sure that material and fill haven't been specified simultaneously.
if ((material[0] != C_NONE) && (fill != C_NONE)) {
if (fill_present && material_present) {
std::stringstream err_msg;
err_msg << "Cell " << id << " has both a material and a fill specified; "
<< "only one can be specified per cell";
fatal_error(err_msg);
}
if (fill_present) {
fill = std::stoi(get_node_value(cell_node, "fill"));
} 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) {
std::vector<std::string> mats
{get_node_array<std::string>(cell_node, "material", true)};
if (mats.size() > 0) {
material.reserve(mats.size());
for (std::string mat : mats) {
if (mat.compare("void") == 0) {
material.push_back(MATERIAL_VOID);
} else {
material.push_back(std::stoi(mat));
}
}
} else {
std::stringstream err_msg;
err_msg << "An empty material element was specified for cell " << id;
fatal_error(err_msg);
}
}
// Read the region specification.
std::string region_spec;
if (check_for_node(cell_node, "region")) {
@ -301,7 +315,7 @@ Cell::distance(Position r, Direction u, int32_t on_surface) const
// Calculate the distance to this surface.
// Note the off-by-one indexing
bool coincident {token == on_surface};
double d {surfaces_c[abs(token)-1]->distance(r, u, coincident)};
double d {global_surfaces[abs(token)-1]->distance(r, u, coincident)};
// Check if this distance is the new minimum.
if (d < min_dist) {
@ -342,7 +356,8 @@ Cell::to_hdf5(hid_t cell_group) const
region_spec << " |";
} else {
// Note the off-by-one indexing
region_spec << " " << copysign(surfaces_c[abs(token)-1]->id, token);
region_spec << " "
<< copysign(global_surfaces[abs(token)-1]->id, token);
}
}
write_string(cell_group, "region", region_spec.str(), false);
@ -365,7 +380,7 @@ Cell::contains_simple(Position r, Direction u, int32_t on_surface) const
return false;
} else {
// Note the off-by-one indexing
bool sense = surfaces_c[abs(token)-1]->sense(r, u);
bool sense = global_surfaces[abs(token)-1]->sense(r, u);
if (sense != (token > 0)) {return false;}
}
}
@ -407,7 +422,7 @@ Cell::contains_complex(Position r, Direction u, int32_t on_surface) const
stack[i_stack] = false;
} else {
// Note the off-by-one indexing
bool sense = surfaces_c[abs(token)-1]->sense(r, u);;
bool sense = global_surfaces[abs(token)-1]->sense(r, u);
stack[i_stack] = (sense == (token > 0));
}
}
@ -429,7 +444,7 @@ Cell::contains_complex(Position r, Direction u, int32_t on_surface) const
//==============================================================================
extern "C" void
read_cells(pugi::xml_node *node)
read_cells(pugi::xml_node* node)
{
// Count the number of cells.
for (pugi::xml_node cell_node: node->children("cell")) {n_cells++;}
@ -437,10 +452,8 @@ read_cells(pugi::xml_node *node)
fatal_error("No cells found in geometry.xml!");
}
// Allocate the vector of Cells.
global_cells.reserve(n_cells);
// Loop over XML cell elements and populate the array.
global_cells.reserve(n_cells);
for (pugi::xml_node cell_node: node->children("cell")) {
global_cells.push_back(new Cell(cell_node));
}
@ -458,6 +471,67 @@ read_cells(pugi::xml_node *node)
global_universes[it->second]->cells.push_back(i);
}
}
global_universes.shrink_to_fit();
}
//==============================================================================
// C-API functions
//==============================================================================
extern "C" int
openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
{
if (index >= 1 && index <= global_cells.size()) {
//TODO: off-by-one
Cell& c {*global_cells[index - 1]};
*type = c.type;
if (c.type == FILL_MATERIAL) {
*indices = c.material.data();
*n = c.material.size();
} else {
*indices = &c.fill;
*n = 1;
}
} else {
strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
}
extern "C" int
openmc_cell_set_fill(int32_t index, int type, int32_t n,
const int32_t* indices)
{
if (index >= 1 && index <= global_cells.size()) {
//TODO: off-by-one
Cell& c {*global_cells[index - 1]};
if (type == FILL_MATERIAL) {
c.type = FILL_MATERIAL;
c.material.clear();
for (int i = 0; i < n; i++) {
int i_mat = indices[i];
if (i_mat == MATERIAL_VOID) {
c.material.push_back(MATERIAL_VOID);
} else if (i_mat >= 1 && i_mat <= global_materials.size()) {
//TODO: off-by-one
c.material.push_back(i_mat - 1);
} else {
strcpy(openmc_err_msg, "Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
c.material.shrink_to_fit();
} else if (type == FILL_UNIVERSE) {
c.type = FILL_UNIVERSE;
} else {
c.type = FILL_LATTICE;
}
} else {
strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
}
//==============================================================================
@ -467,35 +541,39 @@ read_cells(pugi::xml_node *node)
extern "C" {
Cell* cell_pointer(int32_t cell_ind) {return global_cells[cell_ind];}
int32_t cell_id(Cell *c) {return c->id;}
int32_t cell_id(Cell* c) {return c->id;}
void cell_set_id(Cell *c, int32_t id) {c->id = id;}
void cell_set_id(Cell* c, int32_t id) {c->id = id;}
int cell_type(Cell *c) {return c->type;}
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;}
int32_t cell_universe(Cell *c) {return c->universe;}
int32_t cell_fill(Cell* c) {return c->fill;}
void cell_set_universe(Cell *c, int32_t universe) {c->universe = universe;}
int32_t cell_n_instances(Cell* c) {return c->n_instances;}
int32_t cell_fill(Cell *c) {return c->fill;}
int cell_material_size(Cell* c) {return c->material.size();}
int32_t* cell_fill_ptr(Cell *c) {return &c->fill;}
//TODO: off-by-one
int32_t cell_material(Cell* c, int i)
{
int32_t mat = c->material[i-1];
if (mat == MATERIAL_VOID) return MATERIAL_VOID;
return mat + 1;
}
int32_t cell_n_instances(Cell *c) {return c->n_instances;}
bool cell_simple(Cell* c) {return c->simple;}
bool cell_simple(Cell *c) {return c->simple;}
bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface)
bool cell_contains(Cell* c, double xyz[3], double uvw[3], int32_t on_surface)
{
Position r {xyz};
Direction u {uvw};
return c->contains(r, u, on_surface);
}
void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface,
double *min_dist, int32_t *i_surf)
void cell_distance(Cell* c, double xyz[3], double uvw[3], int32_t on_surface,
double* min_dist, int32_t* i_surf)
{
Position r {xyz};
Direction u {uvw};
@ -504,9 +582,9 @@ extern "C" {
*i_surf = out.second;
}
int32_t cell_offset(Cell *c, int map) {return c->offset[map];}
int32_t cell_offset(Cell* c, int map) {return c->offset[map];}
void cell_to_hdf5(Cell *c, hid_t group) {c->to_hdf5(group);}
void cell_to_hdf5(Cell* c, hid_t group) {c->to_hdf5(group);}
void extend_cells_c(int32_t n)
{

View file

@ -172,7 +172,7 @@ contains
p % last_sqrtkT = p % sqrtkT
! Get distributed offset
if (size(c % material) > 1 .or. size(c % sqrtkT) > 1) then
if (c % material_size() > 1 .or. size(c % sqrtkT) > 1) then
! Distributed instances of this cell have different
! materials/temperatures. Determine which instance this is for
! assigning the matching material/temperature.
@ -204,7 +204,7 @@ contains
end if
! Save the material
if (size(c % material) > 1) then
if (c % material_size() > 1) then
p % material = c % material(offset + 1)
else
p % material = c % material(1)

View file

@ -8,6 +8,7 @@
#include "constants.h"
#include "error.h"
#include "lattice.h"
#include "material.h"
namespace openmc {
@ -15,11 +16,11 @@ namespace openmc {
//==============================================================================
void
adjust_indices_c()
adjust_indices()
{
// Adjust material/fill idices.
for (Cell *c : global_cells) {
if (c->material[0] == C_NONE) {
for (Cell* c : global_cells) {
if (c->fill != C_NONE) {
int32_t id = c->fill;
auto search_univ = universe_map.find(id);
auto search_lat = lattice_map.find(id);
@ -36,13 +37,26 @@ adjust_indices_c()
fatal_error(err_msg);
}
} else {
//TODO: materials
c->type = FILL_MATERIAL;
for (auto it = c->material.begin(); it != c->material.end(); it++) {
int32_t mid = *it;
if (mid != MATERIAL_VOID) {
auto search = material_map.find(mid);
if (search != material_map.end()) {
*it = search->second;
} else {
std::stringstream err_msg;
err_msg << "Could not find material " << mid
<< " specified on cell " << c->id;
fatal_error(err_msg);
}
}
}
}
}
// Change cell.universe values from IDs to indices.
for (Cell *c : global_cells) {
for (Cell* c : global_cells) {
auto search = universe_map.find(c->universe);
if (search != universe_map.end()) {
//TODO: Remove this off-by-one indexing.
@ -56,7 +70,7 @@ adjust_indices_c()
}
// Change all lattice universe values from IDs to indices.
for (Lattice *l : lattices_c) {
for (Lattice* l : lattices_c) {
l->adjust_indices();
}
}
@ -68,12 +82,12 @@ find_root_universe()
{
// Find all the universes listed as a cell fill.
std::unordered_set<int32_t> fill_univ_ids;
for (Cell *c : global_cells) {
for (Cell* c : global_cells) {
fill_univ_ids.insert(c->fill);
}
// Find all the universes contained in a lattice.
for (Lattice *lat : lattices_c) {
for (Lattice* lat : lattices_c) {
for (auto it = lat->begin(); it != lat->end(); ++it) {
fill_univ_ids.insert(*it);
}
@ -109,13 +123,13 @@ find_root_universe()
void
allocate_offset_tables(int n_maps)
{
for (Cell *c : global_cells) {
for (Cell* c : global_cells) {
if (c->type != FILL_MATERIAL) {
c->offset.resize(n_maps, C_NONE);
}
}
for (Lattice *lat : lattices_c) {
for (Lattice* lat : lattices_c) {
lat->allocate_offset_table(n_maps);
}
}
@ -126,7 +140,7 @@ void
count_cell_instances(int32_t univ_indx)
{
for (int32_t cell_indx : global_universes[univ_indx]->cells) {
Cell &c = *global_cells[cell_indx];
Cell& c = *global_cells[cell_indx];
++c.n_instances;
if (c.type == FILL_UNIVERSE) {
@ -135,7 +149,7 @@ count_cell_instances(int32_t univ_indx)
} else if (c.type == FILL_LATTICE) {
// This cell contains a lattice. Recurse into the lattice universes.
Lattice &lat = *lattices_c[c.fill];
Lattice& lat = *lattices_c[c.fill];
for (auto it = lat.begin(); it != lat.end(); ++it) {
count_cell_instances(*it);
}
@ -155,14 +169,14 @@ count_universe_instances(int32_t search_univ, int32_t target_univ_id)
int count {0};
for (int32_t cell_indx : global_universes[search_univ]->cells) {
Cell &c = *global_cells[cell_indx];
Cell& c = *global_cells[cell_indx];
if (c.type == FILL_UNIVERSE) {
int32_t next_univ = c.fill;
count += count_universe_instances(next_univ, target_univ_id);
} else if (c.type == FILL_LATTICE) {
Lattice &lat = *lattices_c[c.fill];
Lattice& lat = *lattices_c[c.fill];
for (auto it = lat.begin(); it != lat.end(); ++it) {
int32_t next_univ = *it;
count += count_universe_instances(next_univ, target_univ_id);
@ -178,10 +192,10 @@ count_universe_instances(int32_t search_univ, int32_t target_univ_id)
void
fill_offset_tables(int32_t target_univ_id, int map)
{
for (Universe *univ : global_universes) {
for (Universe* univ : global_universes) {
int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation.
for (int32_t cell_indx : univ->cells) {
Cell &c = *global_cells[cell_indx];
Cell& c = *global_cells[cell_indx];
if (c.type == FILL_UNIVERSE) {
c.offset[map] = offset;
@ -189,7 +203,7 @@ fill_offset_tables(int32_t target_univ_id, int map)
offset += count_universe_instances(search_univ, target_univ_id);
} else if (c.type == FILL_LATTICE) {
Lattice &lat = *lattices_c[c.fill];
Lattice& lat = *lattices_c[c.fill];
offset = lat.fill_offset_table(offset, target_univ_id, map);
}
}
@ -200,7 +214,7 @@ fill_offset_tables(int32_t target_univ_id, int map)
std::string
distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
const Universe &search_univ, int32_t offset)
const Universe& search_univ, int32_t offset)
{
std::stringstream path;
@ -210,7 +224,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
// write to the path and return.
for (int32_t cell_indx : search_univ.cells) {
if ((cell_indx == target_cell) && (offset == target_offset)) {
Cell &c = *global_cells[cell_indx];
Cell& c = *global_cells[cell_indx];
path << "c" << c.id;
return path.str();
}
@ -222,7 +236,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
std::vector<std::int32_t>::const_reverse_iterator cell_it
{search_univ.cells.crbegin()};
for (; cell_it != search_univ.cells.crend(); ++cell_it) {
Cell &c = *global_cells[*cell_it];
Cell& c = *global_cells[*cell_it];
// Material cells don't contain other cells so ignore them.
if (c.type != FILL_MATERIAL) {
@ -230,7 +244,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
if (c.type == FILL_UNIVERSE) {
temp_offset = offset + c.offset[map];
} else {
Lattice &lat = *lattices_c[c.fill];
Lattice& lat = *lattices_c[c.fill];
int32_t indx = lat.universes.size()*map + lat.begin().indx;
temp_offset = offset + lat.offsets[indx];
}
@ -242,7 +256,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
}
// Add the cell to the path string.
Cell &c = *global_cells[*cell_it];
Cell& c = *global_cells[*cell_it];
path << "c" << c.id << "->";
if (c.type == FILL_UNIVERSE) {
@ -253,7 +267,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
return path.str();
} else {
// Recurse into the lattice cell.
Lattice &lat = *lattices_c[c.fill];
Lattice& lat = *lattices_c[c.fill];
path << "l" << lat.id;
for (ReverseLatticeIter it = lat.rbegin(); it != lat.rend(); ++it) {
int32_t indx = lat.universes.size()*map + it.indx;
@ -275,7 +289,7 @@ int
distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset,
int32_t root_univ)
{
Universe &root = *global_universes[root_univ];
Universe& root = *global_universes[root_univ];
std::string path_ {distribcell_path_inner(target_cell, map, target_offset,
root, 0)};
return path_.size() + 1;
@ -285,9 +299,9 @@ distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset,
void
distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset,
int32_t root_univ, char *path)
int32_t root_univ, char* path)
{
Universe &root = *global_universes[root_univ];
Universe& root = *global_universes[root_univ];
std::string path_ {distribcell_path_inner(target_cell, map, target_offset,
root, 0)};
path_.copy(path, path_.size());
@ -302,12 +316,12 @@ maximum_levels(int32_t univ)
int levels_below {0};
for (int32_t cell_indx : global_universes[univ]->cells) {
Cell &c = *global_cells[cell_indx];
Cell& c = *global_cells[cell_indx];
if (c.type == FILL_UNIVERSE) {
int32_t next_univ = c.fill;
levels_below = std::max(levels_below, maximum_levels(next_univ));
} else if (c.type == FILL_LATTICE) {
Lattice &lat = *lattices_c[c.fill];
Lattice& lat = *lattices_c[c.fill];
for (auto it = lat.begin(); it != lat.end(); ++it) {
int32_t next_univ = *it;
levels_below = std::max(levels_below, maximum_levels(next_univ));
@ -324,16 +338,16 @@ maximum_levels(int32_t univ)
void
free_memory_geometry_c()
{
for (Cell *c : global_cells) {delete c;}
for (Cell* c : global_cells) {delete c;}
global_cells.clear();
cell_map.clear();
n_cells = 0;
for (Universe *u : global_universes) {delete u;}
for (Universe* u : global_universes) {delete u;}
global_universes.clear();
universe_map.clear();
for (Lattice *lat : lattices_c) {delete lat;}
for (Lattice* lat : lattices_c) {delete lat;}
lattices_c.clear();
lattice_map.clear();
}

View file

@ -13,7 +13,7 @@ namespace openmc {
//! Replace Universe, Lattice, and Material IDs with indices.
//==============================================================================
extern "C" void adjust_indices_c();
extern "C" void adjust_indices();
//==============================================================================
//! Figure out which Universe is the root universe.

View file

@ -4,7 +4,7 @@ module geometry_header
use algorithm, only: find
use constants, only: HALF, TWO, THREE, INFINITY, K_BOLTZMANN, &
MATERIAL_VOID, NONE
MATERIAL_VOID
use dict_header, only: DictCharInt, DictIntInt
use hdf5_interface, only: HID_T
use material_header, only: Material, materials, material_dict, n_materials
@ -16,11 +16,11 @@ module geometry_header
implicit none
interface
function cell_pointer_c(cell_ind) bind(C, name='cell_pointer') result(ptr)
function cell_pointer(cell_ind) bind(C) result(ptr)
import C_PTR, C_INT32_T
integer(C_INT32_T), intent(in), value :: cell_ind
type(C_PTR) :: ptr
end function cell_pointer_c
end function cell_pointer
function cell_id_c(cell_ptr) bind(C, name='cell_id') result(id)
import C_PTR, C_INT32_T
@ -40,12 +40,6 @@ module geometry_header
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
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
@ -53,25 +47,12 @@ module geometry_header
integer(C_INT32_T) :: universe
end function cell_universe_c
subroutine cell_set_universe_c(cell_ptr, universe) &
bind(C, name='cell_set_universe')
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT32_T), intent(in), value :: universe
end subroutine cell_set_universe_c
function cell_fill_c(cell_ptr) bind(C, name="cell_fill") result(fill)
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT32_T) :: fill
end function cell_fill_c
function cell_fill_ptr(cell_ptr) bind(C) result(fill_ptr)
import C_PTR
type(C_PTR), intent(in), value :: cell_ptr
type(C_PTR) :: fill_ptr
end function cell_fill_ptr
function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') &
result(n_instances)
import C_PTR, C_INT32_T
@ -79,6 +60,21 @@ module geometry_header
integer(C_INT32_T) :: n_instances
end function cell_n_instances_c
function cell_material_size_c(cell_ptr) bind(C, name='cell_material_size') &
result(n)
import C_PTR, C_INT
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT) :: n
end function cell_material_size_c
function cell_material_c(cell_ptr, i) bind(C, name='cell_material') &
result(mat)
import C_PTR, C_INT, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT), intent(in), value :: i
integer(C_INT32_T) :: mat
end function cell_material_c
function cell_simple_c(cell_ptr) bind(C, name='cell_simple') result(simple)
import C_PTR, C_BOOL
type(C_PTR), intent(in), value :: cell_ptr
@ -110,12 +106,11 @@ module geometry_header
integer(HID_T), intent(in), value :: group
end subroutine cell_to_hdf5_c
function lattice_pointer_c(lat_ind) bind(C, name='lattice_pointer') &
result(ptr)
function lattice_pointer(lat_ind) bind(C) result(ptr)
import C_PTR, C_INT32_T
integer(C_INT32_T), intent(in), value :: lat_ind
type(C_PTR) :: ptr
end function lattice_pointer_c
end function lattice_pointer
function lattice_id_c(lat_ptr) bind(C, name='lattice_id') result(id)
import C_PTR, C_INT32_T
@ -254,9 +249,6 @@ module geometry_header
type Cell
type(C_PTR) :: ptr
integer, allocatable :: material(:) ! Material within cell. Multiple
! materials for distribcell
! instances. 0 signifies a universe
integer, allocatable :: region(:) ! Definition of spatial region as
! Boolean expression of half-spaces
integer :: distribcell_index ! Index corresponding to this cell in
@ -275,11 +267,11 @@ 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 :: fill => cell_fill
procedure :: n_instances => cell_n_instances
procedure :: material_size => cell_material_size
procedure :: material => cell_material
procedure :: simple => cell_simple
procedure :: distance => cell_distance
procedure :: offset => cell_offset
@ -390,24 +382,12 @@ contains
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
universe = cell_universe_c(this % ptr)
end function cell_universe
subroutine cell_set_universe(this, universe)
class(Cell), intent(in) :: this
integer(C_INT32_T), intent(in) :: universe
call cell_set_universe_c(this % ptr, universe)
end subroutine cell_set_universe
function cell_fill(this) result(fill)
class(Cell), intent(in) :: this
integer(C_INT32_T) :: fill
@ -420,6 +400,19 @@ contains
n_instances = cell_n_instances_c(this % ptr)
end function cell_n_instances
function cell_material_size(this) result(n)
class(Cell), intent(in) :: this
integer(C_INT) :: n
n = cell_material_size_c(this % ptr)
end function cell_material_size
function cell_material(this, i) result(mat)
class(Cell), intent(in) :: this
integer, intent(in) :: i
integer(C_INT32_T) :: mat
mat = cell_material_c(this % ptr, i)
end function cell_material
function cell_simple(this) result(simple)
class(Cell), intent(in) :: this
logical(C_BOOL) :: simple
@ -469,10 +462,12 @@ contains
if (present(sab_temps)) allocate(sab_temps(n_sab_tables))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
! Skip non-material cells.
if (cells(i) % fill() /= C_NONE) cycle
do j = 1, cells(i) % material_size()
! Skip void materials
if (cells(i) % material(j) == MATERIAL_VOID) cycle
! Get temperature of cell (rounding to nearest integer)
if (size(cells(i) % sqrtkT) > 1) then
@ -571,7 +566,7 @@ contains
! Extend the C++ cells array and get pointers to the C++ objects
call extend_cells_c(n)
do i = n_cells - n, n_cells
cells(i) % ptr = cell_pointer_c(i - 1)
cells(i) % ptr = cell_pointer(i - 1)
end do
err = 0
@ -599,33 +594,6 @@ contains
end function openmc_get_cell_index
function openmc_cell_get_fill(index, type, indices, n) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT), intent(out) :: type
integer(C_INT32_T), intent(out) :: n
type(C_PTR), intent(out) :: indices
integer(C_INT) :: err
err = 0
if (index >= 1 .and. index <= size(cells)) then
associate (c => cells(index))
type = c % type()
select case (type)
case (FILL_MATERIAL)
n = size(c % material)
indices = C_LOC(c % material(1))
case (FILL_UNIVERSE, FILL_LATTICE)
n = 1
indices = cell_fill_ptr(c % ptr)
end select
end associate
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in cells array is out of bounds.")
end if
end function openmc_cell_get_fill
function openmc_cell_get_id(index, id) result(err) bind(C)
! Return the ID of a cell
integer(C_INT32_T), value :: index
@ -642,49 +610,6 @@ contains
end function openmc_cell_get_id
function openmc_cell_set_fill(index, type, n, indices) result(err) bind(C)
! Set the fill for a cell
integer(C_INT32_T), value, intent(in) :: index ! index in cells
integer(C_INT), value, intent(in) :: type
integer(c_INT32_T), value, intent(in) :: n
integer(C_INT32_T), intent(in) :: indices(n)
integer(C_INT) :: err
integer :: i, j
err = 0
if (index >= 1 .and. index <= size(cells)) then
associate (c => cells(index))
select case (type)
case (FILL_MATERIAL)
if (allocated(c % material)) deallocate(c % material)
allocate(c % material(n))
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
c % material(i) = j
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index " // trim(to_str(j)) // " in the &
&materials array is out of bounds.")
end if
end do
case (FILL_UNIVERSE)
call c % set_type(FILL_UNIVERSE)
case (FILL_LATTICE)
call c % set_type(FILL_LATTICE)
end select
end associate
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in cells array is out of bounds.")
end if
end function openmc_cell_set_fill
function openmc_cell_set_id(index, id) result(err) bind(C)
! Set the ID of a cell
integer(C_INT32_T), value, intent(in) :: index

View file

@ -48,8 +48,8 @@ module input_xml
save
interface
subroutine adjust_indices_c() bind(C)
end subroutine adjust_indices_c
subroutine adjust_indices() bind(C)
end subroutine adjust_indices
subroutine allocate_offset_tables(n_maps) bind(C)
import C_INT
@ -87,6 +87,11 @@ module input_xml
type(C_PTR) :: node_ptr
end subroutine read_settings
subroutine read_materials(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine read_materials
function find_root_universe() bind(C) result(root)
import C_INT32_T
integer(C_INT32_T) :: root
@ -1051,7 +1056,7 @@ contains
allocate(surfaces(n_surfaces))
do i = 1, n_surfaces
surfaces(i) % ptr = surface_pointer_c(i - 1);
surfaces(i) % ptr = surface_pointer(i - 1);
if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true.
@ -1095,7 +1100,7 @@ contains
do i = 1, n_cells
c => cells(i)
c % ptr = cell_pointer_c(i - 1)
c % ptr = cell_pointer(i - 1)
! Initialize distribcell instances and distribcell index
c % distribcell_index = NONE
@ -1109,42 +1114,6 @@ contains
// to_str(c % id()))
end if
! Read material
if (check_for_node(node_cell, "material")) then
n_mats = node_word_count(node_cell, "material")
if (n_mats > 0) then
allocate(sarray(n_mats))
call get_node_array(node_cell, "material", sarray)
allocate(c % material(n_mats))
do j = 1, n_mats
select case(trim(to_lower(sarray(j))))
case ('void')
c % material(j) = MATERIAL_VOID
case default
c % material(j) = int(str_to_int(sarray(j)), 4)
! Check for error
if (c % material(j) == ERROR_INT) then
call fatal_error("Invalid material specified on cell " &
// to_str(c % id()))
end if
end select
end do
deallocate(sarray)
else
allocate(c % material(1))
c % material(1) = NONE
end if
else
allocate(c % material(1))
c % material(1) = NONE
end if
! Check for region specification (also under deprecated name surfaces)
if (check_for_node(node_cell, "surfaces")) then
call warning("The use of 'surfaces' is deprecated and will be &
@ -1246,7 +1215,7 @@ contains
n = node_word_count(node_cell, "temperature")
if (n > 0) then
! Make sure this is a "normal" cell.
if (c % material(1) == NONE) call fatal_error("Cell " &
if (c % fill() /= C_NONE) call fatal_error("Cell " &
// trim(to_str(c % id())) // " was specified with a temperature &
&but no material. Temperature specification is only valid for &
&cells filled with a material.")
@ -1309,7 +1278,7 @@ contains
RECT_LATTICES: do i = 1, n_rlats
allocate(RectLattice::lattices(i) % obj)
lat => lattices(i) % obj
lat % ptr = lattice_pointer_c(i - 1)
lat % ptr = lattice_pointer(i - 1)
select type(lat)
type is (RectLattice)
@ -1325,7 +1294,7 @@ contains
HEX_LATTICES: do i = 1, n_hlats
allocate(HexLattice::lattices(n_rlats + i) % obj)
lat => lattices(n_rlats + i) % obj
lat % ptr = lattice_pointer_c(n_rlats + i - 1)
lat % ptr = lattice_pointer(n_rlats + i - 1)
select type (lat)
type is (HexLattice)
@ -1540,10 +1509,12 @@ contains
call doc % load_file(filename)
root = doc % document_element()
call read_materials(root % ptr)
! Get pointer to list of XML <material>
call get_node_list(root, "material", node_mat_list)
! Allocate cells array
! Allocate materials array
n_materials = size(node_mat_list)
allocate(materials(n_materials))
allocate(material_temps(n_materials))
@ -1556,16 +1527,11 @@ contains
do i = 1, n_materials
mat => materials(i)
mat % ptr = material_pointer(i - 1)
! Get pointer to i-th material node
node_mat = node_mat_list(i)
! Copy material id
if (check_for_node(node_mat, "id")) then
call get_node_value(node_mat, "id", mat % id)
else
call fatal_error("Must specify id of material in materials XML file")
end if
! Check if material is depletable
if (check_for_node(node_mat, "depletable")) then
call get_node_value(node_mat, "depletable", temp_str)
@ -1573,12 +1539,6 @@ contains
mat % depletable = .true.
end if
! Check to make sure 'id' hasn't been used
if (material_dict % has(mat % id)) then
call fatal_error("Two or more materials use the same unique ID: " &
// to_str(mat % id))
end if
! Copy material name
if (check_for_node(node_mat, "name")) then
call get_node_value(node_mat, "name", mat % name)
@ -1596,7 +1556,7 @@ contains
node_dens = node_mat % child("density")
else
call fatal_error("Must specify density element in material " &
// trim(to_str(mat % id)))
// trim(to_str(mat % id())))
end if
! Copy units
@ -1626,7 +1586,7 @@ contains
sum_density = .false.
if (val <= ZERO) then
call fatal_error("Need to specify a positive density on material " &
// trim(to_str(mat % id)) // ".")
// trim(to_str(mat % id())) // ".")
end if
! Adjust material density based on specified units
@ -1641,7 +1601,7 @@ contains
mat % density = 1.0e-24_8 * val
case default
call fatal_error("Unkwown units '" // trim(units) &
// "' specified on material " // trim(to_str(mat % id)))
// "' specified on material " // trim(to_str(mat % id())))
end select
end if
@ -1650,7 +1610,7 @@ contains
if (size(node_ele_list) > 0) then
call fatal_error("Unable to add an element to material " &
// trim(to_str(mat % id)) // " since the element option has &
// trim(to_str(mat % id())) // " since the element option has &
&been removed from the xml input. Elements can only be added via &
&the Python API, which will expand elements into their natural &
&nuclides.")
@ -1663,7 +1623,7 @@ contains
if (.not. check_for_node(node_mat, "nuclide") .and. &
.not. check_for_node(node_mat, "macroscopic")) then
call fatal_error("No macroscopic data or nuclides specified on &
&material " // trim(to_str(mat % id)))
&material " // trim(to_str(mat % id())))
end if
! Create list of macroscopic x/s based on those specified, just treat
@ -1677,7 +1637,7 @@ contains
& mode!")
else if (size(node_macro_list) > 1) then
call fatal_error("Only one macroscopic object permitted per material, " &
// trim(to_str(mat % id)))
// trim(to_str(mat % id())))
else if (size(node_macro_list) == 1) then
node_nuc = node_macro_list(1)
@ -1685,7 +1645,7 @@ contains
! Check for empty name on nuclide
if (.not. check_for_node(node_nuc, "name")) then
call fatal_error("No name specified on macroscopic data in material " &
// trim(to_str(mat % id)))
// trim(to_str(mat % id())))
end if
! store nuclide name
@ -1715,7 +1675,7 @@ contains
! Check for empty name on nuclide
if (.not. check_for_node(node_nuc, "name")) then
call fatal_error("No name specified on nuclide in material " &
// trim(to_str(mat % id)))
// trim(to_str(mat % id())))
end if
! store nuclide name
@ -1854,7 +1814,7 @@ contains
if (.not. (all(mat % atom_density >= ZERO) .or. &
all(mat % atom_density <= ZERO))) then
call fatal_error("Cannot mix atom and weight percents in material " &
// to_str(mat % id))
// to_str(mat % id()))
end if
! Determine density if it is a sum value
@ -1935,7 +1895,7 @@ contains
end if
! Add material to dictionary
call material_dict % set(mat % id, i)
call material_dict % set(mat % id(), i)
end do
! Set total number of nuclides and S(a,b) tables
@ -3854,15 +3814,15 @@ contains
do i = 1, n_cells
! Ignore non-normal cells and cells with defined temperature.
if (cells(i) % material(1) == NONE) cycle
if (cells(i) % fill() /= C_NONE) cycle
if (cells(i) % sqrtkT(1) >= ZERO) cycle
! Set the number of temperatures equal to the number of materials.
deallocate(cells(i) % sqrtkT)
allocate(cells(i) % sqrtkT(size(cells(i) % material)))
allocate(cells(i) % sqrtkT(cells(i) % material_size()))
! Check each of the cell materials for temperature data.
do j = 1, size(cells(i) % material)
do j = 1, cells(i) % material_size()
! Arbitrarily set void regions to 0K.
if (cells(i) % material(j) == MATERIAL_VOID) then
cells(i) % sqrtkT(j) = ZERO
@ -3930,45 +3890,6 @@ contains
end subroutine read_multipole_data
!===============================================================================
! ADJUST_INDICES changes the values for 'surfaces' for each cell and the
! material index assigned to each to the indices in the surfaces and material
! array rather than the unique IDs assigned to each surface and material. Also
! assigns boundary conditions to surfaces based on those read into the bc_dict
! dictionary
!===============================================================================
subroutine adjust_indices()
integer :: i ! index for various purposes
integer :: j ! index for various purposes
integer :: id ! user-specified id
call adjust_indices_c()
do i = 1, n_cells
associate (c => cells(i))
! =======================================================================
! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL
if (c % material(1) /= NONE) then
do j = 1, size(c % material)
id = c % material(j)
if (id == MATERIAL_VOID) then
else if (material_dict % has(id)) then
c % material(j) = material_dict % get(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
// " specified on cell " // trim(to_str(c % id())))
end if
end do
end if
end associate
end do
end subroutine adjust_indices
!===============================================================================
! PREPARE_DISTRIBCELL initializes any distribcell filters present and sets the
! offsets for distribcells
@ -3992,7 +3913,7 @@ contains
! Find all cells with multiple (distributed) materials or temperatures.
do i = 1, n_cells
if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then
if (cells(i) % material_size() > 1 .or. size(cells(i) % sqrtkT) > 1) then
call cell_list % add(i)
end if
end do
@ -4001,10 +3922,10 @@ contains
! number of respective cell instances.
do i = 1, n_cells
associate (c => cells(i))
if (size(c % material) > 1) then
if (size(c % material) /= c % n_instances()) then
if (c % material_size() > 1) then
if (c % material_size() /= c % n_instances()) then
call fatal_error("Cell " // trim(to_str(c % id())) // " was &
&specified with " // trim(to_str(size(c % material))) &
&specified with " // trim(to_str(c % material_size())) &
// " materials but has " // trim(to_str(c % n_instances())) &
// " distributed instances. The number of materials must &
&equal one or the number of instances.")

View file

@ -29,7 +29,7 @@ std::unordered_map<int32_t, int32_t> lattice_map;
Lattice::Lattice(pugi::xml_node lat_node)
{
if (check_for_node(lat_node, "id")) {
id = stoi(get_node_value(lat_node, "id"));
id = std::stoi(get_node_value(lat_node, "id"));
} else {
fatal_error("Must specify id of lattice in geometry XML file.");
}
@ -39,7 +39,7 @@ Lattice::Lattice(pugi::xml_node lat_node)
}
if (check_for_node(lat_node, "outer")) {
outer = stoi(get_node_value(lat_node, "outer"));
outer = std::stoi(get_node_value(lat_node, "outer"));
}
}
@ -141,14 +141,14 @@ RectLattice::RectLattice(pugi::xml_node lat_node)
std::string dimension_str {get_node_value(lat_node, "dimension")};
std::vector<std::string> dimension_words {split(dimension_str)};
if (dimension_words.size() == 2) {
n_cells[0] = stoi(dimension_words[0]);
n_cells[1] = stoi(dimension_words[1]);
n_cells[0] = std::stoi(dimension_words[0]);
n_cells[1] = std::stoi(dimension_words[1]);
n_cells[2] = 1;
is_3d = false;
} else if (dimension_words.size() == 3) {
n_cells[0] = stoi(dimension_words[0]);
n_cells[1] = stoi(dimension_words[1]);
n_cells[2] = stoi(dimension_words[2]);
n_cells[0] = std::stoi(dimension_words[0]);
n_cells[1] = std::stoi(dimension_words[1]);
n_cells[2] = std::stoi(dimension_words[2]);
is_3d = true;
} else {
fatal_error("Rectangular lattice must be two or three dimensions.");
@ -195,7 +195,7 @@ RectLattice::RectLattice(pugi::xml_node lat_node)
for (int ix = 0; ix < nx; ix++) {
int indx1 = nx*ny*iz + nx*(ny-iy-1) + ix;
int indx2 = nx*ny*iz + nx*iy + ix;
universes[indx1] = stoi(univ_words[indx2]);
universes[indx1] = std::stoi(univ_words[indx2]);
}
}
}
@ -400,9 +400,9 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
: Lattice {lat_node}
{
// Read the number of lattice cells in each dimension.
n_rings = stoi(get_node_value(lat_node, "n_rings"));
n_rings = std::stoi(get_node_value(lat_node, "n_rings"));
if (check_for_node(lat_node, "n_axial")) {
n_axial = stoi(get_node_value(lat_node, "n_axial"));
n_axial = std::stoi(get_node_value(lat_node, "n_axial"));
is_3d = true;
} else {
n_axial = 1;
@ -476,7 +476,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
int indx = (2*n_rings-1)*(2*n_rings-1) * m
+ (2*n_rings-1) * (i_a+n_rings-1)
+ (i_x+n_rings-1);
universes[indx] = stoi(univ_words[input_index]);
universes[indx] = std::stoi(univ_words[input_index]);
input_index++;
// Walk the index to the right neighbor (which is not adjacent).
i_x += 2;
@ -505,7 +505,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
int indx = (2*n_rings-1)*(2*n_rings-1) * m
+ (2*n_rings-1) * (i_a+n_rings-1)
+ (i_x+n_rings-1);
universes[indx] = stoi(univ_words[input_index]);
universes[indx] = std::stoi(univ_words[input_index]);
input_index++;
// Walk the index to the right neighbor (which is not adjacent).
i_x += 2;
@ -528,7 +528,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
int indx = (2*n_rings-1)*(2*n_rings-1) * m
+ (2*n_rings-1) * (i_a+n_rings-1)
+ (i_x+n_rings-1);
universes[indx] = stoi(univ_words[input_index]);
universes[indx] = std::stoi(univ_words[input_index]);
input_index++;
// Walk the index to the right neighbor (which is not adjacent).
i_x += 2;

91
src/material.cpp Normal file
View file

@ -0,0 +1,91 @@
#include "material.h"
#include <string>
#include <sstream>
#include "error.h"
#include "xml_interface.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
std::vector<Material*> global_materials;
std::unordered_map<int32_t, int32_t> material_map;
//==============================================================================
// Material implementation
//==============================================================================
Material::Material(pugi::xml_node material_node)
{
if (check_for_node(material_node, "id")) {
id = std::stoi(get_node_value(material_node, "id"));
} else {
fatal_error("Must specify id of material in materials XML file.");
}
}
//==============================================================================
// Non-method functions
//==============================================================================
extern "C" void
read_materials(pugi::xml_node* node)
{
// Loop over XML material elements and populate the array.
for (pugi::xml_node material_node : node->children("material")) {
global_materials.push_back(new Material(material_node));
}
global_materials.shrink_to_fit();
// Populate the material map.
for (int i = 0; i < global_materials.size(); i++) {
int32_t mid = global_materials[i]->id;
auto search = material_map.find(mid);
if (search == material_map.end()) {
material_map[mid] = i;
} else {
std::stringstream err_msg;
err_msg << "Two or more materials use the same unique ID: " << mid;
fatal_error(err_msg);
}
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
Material* material_pointer(int32_t indx) {return global_materials[indx];}
int32_t material_id(Material* mat) {return mat->id;}
void material_set_id(Material* mat, int32_t id, int32_t index)
{
mat->id = id;
//TODO: off-by-one
material_map[id] = index - 1;
}
void extend_materials_c(int32_t n)
{
global_materials.reserve(global_materials.size() + n);
for (int32_t i = 0; i < n; i++) {
global_materials.push_back(new Material());
}
}
void free_memory_material_c()
{
for (Material *mat : global_materials) {delete mat;}
global_materials.clear();
material_map.clear();
}
}
} // namespace openmc

35
src/material.h Normal file
View file

@ -0,0 +1,35 @@
#ifndef OPENMC_MATERIAL_H
#define OPENMC_MATERIAL_H
#include <unordered_map>
#include <vector>
#include "pugixml.hpp"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
class Material;
extern std::vector<Material*> global_materials;
extern std::unordered_map<int32_t, int32_t> material_map;
//==============================================================================
//! A substance with constituent nuclides and thermal scattering data
//==============================================================================
class Material
{
public:
int32_t id; //!< Unique ID
Material() {};
explicit Material(pugi::xml_node material_node);
};
} // namespace openmc
#endif // OPENMC_MATERIAL_H

View file

@ -27,13 +27,41 @@ module material_header
public :: openmc_material_set_density
public :: openmc_material_set_densities
public :: openmc_material_set_id
public :: material_pointer
interface
function material_pointer(mat_ind) bind(C) result(ptr)
import C_PTR, C_INT32_T
integer(C_INT32_T), intent(in), value :: mat_ind
type(C_PTR) :: ptr
end function material_pointer
function material_id_c(mat_ptr) bind(C, name='material_id') result(id)
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: mat_ptr
integer(C_INT32_T) :: id
end function material_id_c
subroutine material_set_id_c(mat_ptr, id, index) &
bind(C, name='material_set_id')
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: mat_ptr
integer(C_INT32_T), intent(in), value :: id
integer(C_INT32_T), intent(in), value :: index
end subroutine material_set_id_c
subroutine extend_materials_c(n) bind(C)
import C_INT32_t
integer(C_INT32_T), intent(in), value :: n
end subroutine extend_materials_c
end interface
!===============================================================================
! MATERIAL describes a material by its constituent nuclides
!===============================================================================
type, public :: Material
integer :: id ! unique identifier
type(C_PTR) :: ptr
character(len=104) :: name = "" ! User-defined name
integer :: n_nuclides = 0 ! number of nuclides
integer, allocatable :: nuclide(:) ! index in nuclides array
@ -67,6 +95,8 @@ module material_header
logical, allocatable :: p0(:)
contains
procedure :: id => material_id
procedure :: set_id => material_set_id
procedure :: set_density => material_set_density
procedure :: init_nuclide_index => material_init_nuclide_index
procedure :: assign_sab_tables => material_assign_sab_tables
@ -88,6 +118,19 @@ contains
! MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm.
!===============================================================================
function material_id(this) result(id)
class(Material), intent(in) :: this
integer(C_INT32_T) :: id
id = material_id_c(this % ptr)
end function material_id
subroutine material_set_id(this, id, index)
class(Material), intent(in) :: this
integer(C_INT32_T), intent(in) :: id
integer(C_INT32_T), intent(in) :: index
call material_set_id_c(this % ptr, id, index)
end subroutine material_set_id
function material_set_density(this, density) result(err)
class(Material), intent(inout) :: this
real(8), intent(in) :: density
@ -185,7 +228,7 @@ contains
if (.not. found) then
call fatal_error("S(a,b) table " // trim(this % &
sab_names(k)) // " did not match any nuclide on material " &
// trim(to_str(this % id)))
// trim(to_str(this % id())))
end if
end do ASSIGN_SAB
@ -195,7 +238,7 @@ contains
if (i_sab_nuclides % data(j) == i_sab_nuclides % data(k)) then
call fatal_error(trim( &
nuclides(this % nuclide(i_sab_nuclides % data(j))) % name) &
// " in material " // trim(to_str(this % id)) // " was found &
// " in material " // trim(to_str(this % id())) // " was found &
&in multiple S(a,b) tables. Each nuclide can only appear in &
&one S(a,b) table per material.")
end if
@ -444,6 +487,11 @@ contains
!===============================================================================
subroutine free_memory_material()
interface
subroutine free_memory_material_c() bind(C)
end subroutine free_memory_material_c
end interface
call free_memory_material_c()
n_materials = 0
if (allocated(materials)) deallocate(materials)
call material_dict % clear()
@ -460,6 +508,7 @@ contains
integer(C_INT32_T), optional, intent(out) :: index_end
integer(C_INT) :: err
integer :: i
type(Material), allocatable :: temp(:) ! temporary materials array
if (n_materials == 0) then
@ -481,6 +530,12 @@ contains
if (present(index_end)) index_end = n_materials + n
n_materials = n_materials + n
! Extend the C++ materials array and get pointers to the C++ objects
call extend_materials_c(n)
do i = n_materials - n, n_materials
materials(i) % ptr = material_pointer(i - 1)
end do
err = 0
end function openmc_extend_materials
@ -606,7 +661,7 @@ contains
integer(C_INT) :: err
if (index >= 1 .and. index <= size(materials)) then
id = materials(index) % id
id = materials(index) % id()
err = 0
else
err = E_OUT_OF_BOUNDS
@ -622,7 +677,7 @@ contains
integer(C_INT) :: err
if (index >= 1 .and. index <= n_materials) then
materials(index) % id = id
call materials(index) % set_id(id, index)
call material_dict % set(id, index)
err = 0
else

View file

@ -151,11 +151,13 @@ contains
allocate(kTs(size(materials)))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip non-material cells
if (cells(i) % fill() /= C_NONE) cycle
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
do j = 1, cells(i) % material_size()
! Skip void materials
if (cells(i) % material(j) == MATERIAL_VOID) cycle
! Get temperature of cell (rounding to nearest integer)
if (size(cells(i) % sqrtkT) > 1) then

View file

@ -95,7 +95,7 @@ contains
id = -1
else
rgb = pl % colors(p % material) % rgb
id = materials(p % material) % id
id = materials(p % material) % id()
end if
end associate
else if (pl % color_by == PLOT_COLOR_CELLS) then

View file

@ -198,20 +198,20 @@ contains
case (FILL_MATERIAL)
call write_dataset(cell_group, "fill_type", "material")
if (size(c % material) == 1) then
if (c % material_size() == 1) then
if (c % material(1) == MATERIAL_VOID) then
call write_dataset(cell_group, "material", MATERIAL_VOID)
else
call write_dataset(cell_group, "material", &
materials(c % material(1)) % id)
materials(c % material(1)) % id())
end if
else
allocate(cell_materials(size(c % material)))
do j = 1, size(c % material)
allocate(cell_materials(c % material_size()))
do j = 1, c % material_size()
if (c % material(j) == MATERIAL_VOID) then
cell_materials(j) = MATERIAL_VOID
else
cell_materials(j) = materials(c % material(j)) % id
cell_materials(j) = materials(c % material(j)) % id()
end if
end do
call write_dataset(cell_group, "material", cell_materials)
@ -333,7 +333,7 @@ contains
do i = 1, n_materials
m => materials(i)
material_group = create_group(materials_group, "material " // &
trim(to_str(m%id)))
trim(to_str(m%id())))
if (m % depletable) then
call write_attribute(material_group, "depletable", 1)

View file

@ -27,7 +27,7 @@ extern "C" const int BC_PERIODIC {3};
int32_t n_surfaces;
Surface **surfaces_c;
std::vector<Surface*> global_surfaces;
std::map<int, int> surface_map;
@ -141,13 +141,13 @@ void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1, double &c2,
Surface::Surface(pugi::xml_node surf_node)
{
if (check_for_node(surf_node, "id")) {
id = stoi(get_node_value(surf_node, "id"));
id = std::stoi(get_node_value(surf_node, "id"));
} else {
fatal_error("Must specify id of surface in geometry XML file.");
}
if (check_for_node(surf_node, "name")) {
name = get_node_value(surf_node, "name");
name = get_node_value(surf_node, "name", false);
}
if (check_for_node(surf_node, "boundary")) {
@ -247,7 +247,7 @@ PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node)
: Surface {surf_node}
{
if (check_for_node(surf_node, "periodic_surface_id")) {
i_periodic = stoi(get_node_value(surf_node, "periodic_surface_id"));
i_periodic = std::stoi(get_node_value(surf_node, "periodic_surface_id"));
}
}
@ -298,8 +298,8 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfaceXPlane::periodic_translate(const PeriodicSurface *other, Position& r,
Direction& u) const
bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other,
Position& r, Direction& u) const
{
Direction other_n = other->normal(r);
if (other_n.x == 1 and other_n.y == 0 and other_n.z == 0) {
@ -359,8 +359,8 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, Position& r,
Direction& u) const
bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other,
Position& r, Direction& u) const
{
Direction other_n = other->normal(r);
if (other_n.x == 0 and other_n.y == 1 and other_n.z == 0) {
@ -421,8 +421,8 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfaceZPlane::periodic_translate(const PeriodicSurface *other, Position& r,
Direction& u) const
bool SurfaceZPlane::periodic_translate(const PeriodicSurface* other,
Position& r, Direction& u) const
{
// Assume the other plane is aligned along z. Just change the z coord.
r.z = z0;
@ -478,7 +478,7 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfacePlane::periodic_translate(const PeriodicSurface *other, Position& r,
bool SurfacePlane::periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const
{
// This function assumes the other plane shares this plane's normal direction.
@ -1023,7 +1023,7 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
//==============================================================================
extern "C" void
read_surfaces(pugi::xml_node *node)
read_surfaces(pugi::xml_node* node)
{
// Count the number of surfaces.
for (pugi::xml_node surf_node: node->children("surface")) {n_surfaces++;}
@ -1031,10 +1031,8 @@ read_surfaces(pugi::xml_node *node)
fatal_error("No surfaces found in geometry.xml!");
}
// Allocate the array of Surface pointers.
surfaces_c = new Surface* [n_surfaces];
// Loop over XML surface elements and populate the array.
global_surfaces.reserve(n_surfaces);
{
pugi::xml_node surf_node;
int i_surf;
@ -1043,40 +1041,40 @@ read_surfaces(pugi::xml_node *node)
std::string surf_type = get_node_value(surf_node, "type", true, true);
if (surf_type == "x-plane") {
surfaces_c[i_surf] = new SurfaceXPlane(surf_node);
global_surfaces.push_back(new SurfaceXPlane(surf_node));
} else if (surf_type == "y-plane") {
surfaces_c[i_surf] = new SurfaceYPlane(surf_node);
global_surfaces.push_back(new SurfaceYPlane(surf_node));
} else if (surf_type == "z-plane") {
surfaces_c[i_surf] = new SurfaceZPlane(surf_node);
global_surfaces.push_back(new SurfaceZPlane(surf_node));
} else if (surf_type == "plane") {
surfaces_c[i_surf] = new SurfacePlane(surf_node);
global_surfaces.push_back(new SurfacePlane(surf_node));
} else if (surf_type == "x-cylinder") {
surfaces_c[i_surf] = new SurfaceXCylinder(surf_node);
global_surfaces.push_back(new SurfaceXCylinder(surf_node));
} else if (surf_type == "y-cylinder") {
surfaces_c[i_surf] = new SurfaceYCylinder(surf_node);
global_surfaces.push_back(new SurfaceYCylinder(surf_node));
} else if (surf_type == "z-cylinder") {
surfaces_c[i_surf] = new SurfaceZCylinder(surf_node);
global_surfaces.push_back(new SurfaceZCylinder(surf_node));
} else if (surf_type == "sphere") {
surfaces_c[i_surf] = new SurfaceSphere(surf_node);
global_surfaces.push_back(new SurfaceSphere(surf_node));
} else if (surf_type == "x-cone") {
surfaces_c[i_surf] = new SurfaceXCone(surf_node);
global_surfaces.push_back(new SurfaceXCone(surf_node));
} else if (surf_type == "y-cone") {
surfaces_c[i_surf] = new SurfaceYCone(surf_node);
global_surfaces.push_back(new SurfaceYCone(surf_node));
} else if (surf_type == "z-cone") {
surfaces_c[i_surf] = new SurfaceZCone(surf_node);
global_surfaces.push_back(new SurfaceZCone(surf_node));
} else if (surf_type == "quadric") {
surfaces_c[i_surf] = new SurfaceQuadric(surf_node);
global_surfaces.push_back(new SurfaceQuadric(surf_node));
} else {
std::stringstream err_msg;
@ -1088,7 +1086,7 @@ read_surfaces(pugi::xml_node *node)
// Fill the surface map.
for (int i_surf = 0; i_surf < n_surfaces; i_surf++) {
int id = surfaces_c[i_surf]->id;
int id = global_surfaces[i_surf]->id;
auto in_map = surface_map.find(id);
if (in_map == surface_map.end()) {
surface_map[id] = i_surf;
@ -1104,10 +1102,10 @@ read_surfaces(pugi::xml_node *node)
zmin {INFTY}, zmax {-INFTY};
int i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax;
for (int i_surf = 0; i_surf < n_surfaces; i_surf++) {
if (surfaces_c[i_surf]->bc == BC_PERIODIC) {
if (global_surfaces[i_surf]->bc == BC_PERIODIC) {
// Downcast to the PeriodicSurface type.
Surface *surf_base = surfaces_c[i_surf];
PeriodicSurface *surf = dynamic_cast<PeriodicSurface *>(surf_base);
Surface* surf_base = global_surfaces[i_surf];
PeriodicSurface* surf = dynamic_cast<PeriodicSurface*>(surf_base);
// Make sure this surface inherits from PeriodicSurface.
if (!surf) {
@ -1149,14 +1147,14 @@ read_surfaces(pugi::xml_node *node)
// Set i_periodic for periodic BC surfaces.
for (int i_surf = 0; i_surf < n_surfaces; i_surf++) {
if (surfaces_c[i_surf]->bc == BC_PERIODIC) {
if (global_surfaces[i_surf]->bc == BC_PERIODIC) {
// Downcast to the PeriodicSurface type.
Surface *surf_base = surfaces_c[i_surf];
PeriodicSurface *surf = dynamic_cast<PeriodicSurface *>(surf_base);
Surface* surf_base = global_surfaces[i_surf];
PeriodicSurface* surf = dynamic_cast<PeriodicSurface*>(surf_base);
// Also try downcasting to the SurfacePlane type (which must be handled
// differently).
SurfacePlane *surf_p = dynamic_cast<SurfacePlane *>(surf);
SurfacePlane* surf_p = dynamic_cast<SurfacePlane*>(surf);
if (!surf_p) {
// This is not a SurfacePlane.
@ -1198,7 +1196,7 @@ read_surfaces(pugi::xml_node *node)
}
// Make sure the opposite surface is also periodic.
if (surfaces_c[surf->i_periodic]->bc != BC_PERIODIC) {
if (global_surfaces[surf->i_periodic]->bc != BC_PERIODIC) {
std::stringstream err_msg;
err_msg << "Could not find matching surface for periodic boundary "
"condition on surface " << surf->id;
@ -1213,13 +1211,13 @@ read_surfaces(pugi::xml_node *node)
//==============================================================================
extern "C" {
Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];}
Surface* surface_pointer(int surf_ind) {return global_surfaces[surf_ind];}
int surface_id(Surface *surf) {return surf->id;}
int surface_id(Surface* surf) {return surf->id;}
int surface_bc(Surface *surf) {return surf->bc;}
int surface_bc(Surface* surf) {return surf->bc;}
void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
void surface_reflect(Surface* surf, double xyz[3], double uvw[3])
{
Position r {xyz};
Direction u {uvw};
@ -1230,7 +1228,7 @@ extern "C" {
uvw[2] = u.z;
}
void surface_normal(Surface *surf, double xyz[3], double uvw[3])
void surface_normal(Surface* surf, double xyz[3], double uvw[3])
{
Position r {xyz};
Direction u = surf->normal(r);
@ -1239,12 +1237,12 @@ extern "C" {
uvw[2] = u.z;
}
void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);}
void surface_to_hdf5(Surface* surf, hid_t group) {surf->to_hdf5(group);}
int surface_i_periodic(PeriodicSurface *surf) {return surf->i_periodic;}
int surface_i_periodic(PeriodicSurface* surf) {return surf->i_periodic;}
bool
surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3],
surface_periodic(PeriodicSurface* surf, PeriodicSurface* other, double xyz[3],
double uvw[3])
{
Position r {xyz};
@ -1264,9 +1262,8 @@ extern "C" {
void free_memory_surfaces_c()
{
for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];}
delete surfaces_c;
surfaces_c = nullptr;
for (Surface* surf : global_surfaces) {delete surf;}
global_surfaces.clear();
n_surfaces = 0;
surface_map.clear();
}

View file

@ -4,6 +4,7 @@
#include <map>
#include <limits> // For numeric_limits
#include <string>
#include <vector>
#include "hdf5.h"
#include "pugixml.hpp"
@ -31,7 +32,7 @@ extern "C" const int BC_PERIODIC;
extern "C" int32_t n_surfaces;
class Surface;
extern Surface **surfaces_c;
extern std::vector<Surface*> global_surfaces;
extern std::map<int, int> surface_map;
@ -131,7 +132,7 @@ public:
//! periodicity.
//! \return true if this surface and its partner make a rotationally-periodic
//! boundary condition.
virtual bool periodic_translate(const PeriodicSurface *other, Position& r,
virtual bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const = 0;
//! Get the bounding box for this surface.
@ -153,8 +154,8 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
};
@ -173,8 +174,8 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
};
@ -193,8 +194,8 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
};
@ -213,8 +214,8 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
};
@ -368,16 +369,16 @@ public:
extern "C" {
Surface* surface_pointer(int surf_ind);
int surface_id(Surface *surf);
int surface_bc(Surface *surf);
bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
double surface_distance(Surface *surf, double xyz[3], double uvw[3],
int surface_id(Surface* surf);
int surface_bc(Surface* surf);
bool surface_sense(Surface* surf, double xyz[3], double uvw[3]);
void surface_reflect(Surface* surf, double xyz[3], double uvw[3]);
double surface_distance(Surface* surf, double xyz[3], double uvw[3],
bool coincident);
void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
void surface_to_hdf5(Surface *surf, hid_t group);
int surface_i_periodic(PeriodicSurface *surf);
bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
void surface_normal(Surface* surf, double xyz[3], double uvw[3]);
void surface_to_hdf5(Surface* surf, hid_t group);
int surface_i_periodic(PeriodicSurface* surf);
bool surface_periodic(PeriodicSurface* surf, PeriodicSurface* other,
double xyz[3], double uvw[3]);
void free_memory_surfaces_c();
}

View file

@ -8,13 +8,12 @@ module surface_header
implicit none
interface
pure function surface_pointer_c(surf_ind) &
bind(C, name='surface_pointer') result(ptr)
pure function surface_pointer(surf_ind) bind(C) result(ptr)
use ISO_C_BINDING
implicit none
integer(C_INT), intent(in), value :: surf_ind
type(C_PTR) :: ptr
end function surface_pointer_c
end function surface_pointer
pure function surface_id_c(surf_ptr) bind(C, name='surface_id') result(id)
use ISO_C_BINDING

View file

@ -3043,7 +3043,7 @@ contains
case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, &
SCORE_NU_FISSION)
if (materials(p % material) % id == deriv % diff_material) then
if (materials(p % material) % id() == deriv % diff_material) then
score = score * (flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
@ -3064,7 +3064,7 @@ contains
case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, &
SCORE_NU_FISSION)
if (materials(p % material) % id == deriv % diff_material) then
if (materials(p % material) % id() == deriv % diff_material) then
score = score * (flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
@ -3106,7 +3106,7 @@ contains
case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, &
SCORE_NU_FISSION)
if (materials(p % material) % id == deriv % diff_material &
if (materials(p % material) % id() == deriv % diff_material &
.and. p % event_nuclide == deriv % diff_nuclide) then
associate(mat => materials(p % material))
! Search for the index of the perturbed nuclide.
@ -3128,7 +3128,7 @@ contains
case (ESTIMATOR_COLLISION)
scoring_diff_nuclide = &
(materials(p % material) % id == deriv % diff_material) &
(materials(p % material) % id() == deriv % diff_material) &
.and. (i_nuclide == deriv % diff_nuclide)
select case (score_bin)
@ -3138,7 +3138,7 @@ contains
case (SCORE_TOTAL)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material .and. &
materials(p % material) % id() == deriv % diff_material .and. &
material_xs % total /= ZERO) then
score = score * (flux_deriv &
+ micro_xs(deriv % diff_nuclide) % total &
@ -3152,7 +3152,7 @@ contains
case (SCORE_SCATTER)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material .and. &
materials(p % material) % id() == deriv % diff_material .and. &
material_xs % total - material_xs % absorption /= ZERO) then
score = score * (flux_deriv &
+ (micro_xs(deriv % diff_nuclide) % total &
@ -3168,7 +3168,7 @@ contains
case (SCORE_ABSORPTION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material .and. &
materials(p % material) % id() == deriv % diff_material .and. &
material_xs % absorption /= ZERO) then
score = score * (flux_deriv &
+ micro_xs(deriv % diff_nuclide) % absorption &
@ -3182,7 +3182,7 @@ contains
case (SCORE_FISSION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material .and. &
materials(p % material) % id() == deriv % diff_material .and. &
material_xs % fission /= ZERO) then
score = score * (flux_deriv &
+ micro_xs(deriv % diff_nuclide) % fission &
@ -3196,7 +3196,7 @@ contains
case (SCORE_NU_FISSION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material .and. &
materials(p % material) % id() == deriv % diff_material .and. &
material_xs % nu_fission /= ZERO) then
score = score * (flux_deriv &
+ micro_xs(deriv % diff_nuclide) % nu_fission &
@ -3242,7 +3242,7 @@ contains
score = score * flux_deriv
case (SCORE_TOTAL)
if (materials(p % material) % id == deriv % diff_material .and. &
if (materials(p % material) % id() == deriv % diff_material .and. &
micro_xs(p % event_nuclide) % total > ZERO) then
associate(mat => materials(p % material))
! Search for the index of the perturbed nuclide.
@ -3267,7 +3267,7 @@ contains
end if
case (SCORE_SCATTER)
if (materials(p % material) % id == deriv % diff_material .and. &
if (materials(p % material) % id() == deriv % diff_material .and. &
(micro_xs(p % event_nuclide) % total &
- micro_xs(p % event_nuclide) % absorption) > ZERO) then
associate(mat => materials(p % material))
@ -3295,7 +3295,7 @@ contains
end if
case (SCORE_ABSORPTION)
if (materials(p % material) % id == deriv % diff_material .and. &
if (materials(p % material) % id() == deriv % diff_material .and. &
micro_xs(p % event_nuclide) % absorption > ZERO) then
associate(mat => materials(p % material))
! Search for the index of the perturbed nuclide.
@ -3320,7 +3320,7 @@ contains
end if
case (SCORE_FISSION)
if (materials(p % material) % id == deriv % diff_material .and. &
if (materials(p % material) % id() == deriv % diff_material .and. &
micro_xs(p % event_nuclide) % fission > ZERO) then
associate(mat => materials(p % material))
! Search for the index of the perturbed nuclide.
@ -3345,7 +3345,7 @@ contains
end if
case (SCORE_NU_FISSION)
if (materials(p % material) % id == deriv % diff_material .and. &
if (materials(p % material) % id() == deriv % diff_material .and. &
micro_xs(p % event_nuclide) % nu_fission > ZERO) then
associate(mat => materials(p % material))
! Search for the index of the perturbed nuclide.
@ -3385,7 +3385,7 @@ contains
case (SCORE_TOTAL)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material .and. &
materials(p % material) % id() == deriv % diff_material .and. &
material_xs % total > ZERO) then
cum_dsig = ZERO
associate(mat => materials(p % material))
@ -3404,7 +3404,7 @@ contains
end associate
score = score * (flux_deriv &
+ cum_dsig / material_xs % total)
else if (materials(p % material) % id == deriv % diff_material &
else if (materials(p % material) % id() == deriv % diff_material &
.and. material_xs % total > ZERO) then
dsig_t = ZERO
associate (nuc => nuclides(i_nuclide))
@ -3423,7 +3423,7 @@ contains
case (SCORE_SCATTER)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material .and. &
materials(p % material) % id() == deriv % diff_material .and. &
(material_xs % total - material_xs % absorption) > ZERO) then
cum_dsig = ZERO
associate(mat => materials(p % material))
@ -3444,7 +3444,7 @@ contains
end associate
score = score * (flux_deriv + cum_dsig &
/ (material_xs % total - material_xs % absorption))
else if ( materials(p % material) % id == deriv % diff_material &
else if ( materials(p % material) % id() == deriv % diff_material &
.and. (material_xs % total - material_xs % absorption) > ZERO)&
then
dsig_t = ZERO
@ -3466,7 +3466,7 @@ contains
case (SCORE_ABSORPTION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material .and. &
materials(p % material) % id() == deriv % diff_material .and. &
material_xs % absorption > ZERO) then
cum_dsig = ZERO
associate(mat => materials(p % material))
@ -3485,7 +3485,7 @@ contains
end associate
score = score * (flux_deriv &
+ cum_dsig / material_xs % absorption)
else if (materials(p % material) % id == deriv % diff_material &
else if (materials(p % material) % id() == deriv % diff_material &
.and. material_xs % absorption > ZERO) then
dsig_a = ZERO
associate (nuc => nuclides(i_nuclide))
@ -3504,7 +3504,7 @@ contains
case (SCORE_FISSION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material .and. &
materials(p % material) % id() == deriv % diff_material .and. &
material_xs % fission > ZERO) then
cum_dsig = ZERO
associate(mat => materials(p % material))
@ -3523,7 +3523,7 @@ contains
end associate
score = score * (flux_deriv &
+ cum_dsig / material_xs % fission)
else if (materials(p % material) % id == deriv % diff_material &
else if (materials(p % material) % id() == deriv % diff_material &
.and. material_xs % fission > ZERO) then
dsig_f = ZERO
associate (nuc => nuclides(i_nuclide))
@ -3542,7 +3542,7 @@ contains
case (SCORE_NU_FISSION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material .and. &
materials(p % material) % id() == deriv % diff_material .and. &
material_xs % nu_fission > ZERO) then
cum_dsig = ZERO
associate(mat => materials(p % material))
@ -3563,7 +3563,7 @@ contains
end associate
score = score * (flux_deriv &
+ cum_dsig / material_xs % nu_fission)
else if (materials(p % material) % id == deriv % diff_material &
else if (materials(p % material) % id() == deriv % diff_material &
.and. material_xs % nu_fission > ZERO) then
dsig_f = ZERO
associate (nuc => nuclides(i_nuclide))
@ -3614,7 +3614,7 @@ contains
case (DIFF_DENSITY)
associate (mat => materials(p % material))
if (mat % id == deriv % diff_material) then
if (mat % id() == deriv % diff_material) then
! phi is proportional to e^(-Sigma_tot * dist)
! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
@ -3625,7 +3625,7 @@ contains
case (DIFF_NUCLIDE_DENSITY)
associate (mat => materials(p % material))
if (mat % id == deriv % diff_material) then
if (mat % id() == deriv % diff_material) then
! phi is proportional to e^(-Sigma_tot * dist)
! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist
! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist
@ -3636,7 +3636,7 @@ contains
case (DIFF_TEMPERATURE)
associate (mat => materials(p % material))
if (mat % id == deriv % diff_material) then
if (mat % id() == deriv % diff_material) then
do l=1, mat % n_nuclides
associate (nuc => nuclides(mat % nuclide(l)))
if (nuc % mp_present .and. &
@ -3690,7 +3690,7 @@ contains
case (DIFF_DENSITY)
associate (mat => materials(p % material))
if (mat % id == deriv % diff_material) then
if (mat % id() == deriv % diff_material) then
! phi is proportional to Sigma_s
! (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s
! (1 / phi) * (d_phi / d_rho) = 1 / rho
@ -3701,7 +3701,7 @@ contains
case (DIFF_NUCLIDE_DENSITY)
associate (mat => materials(p % material))
if (mat % id == deriv % diff_material &
if (mat % id() == deriv % diff_material &
.and. p % event_nuclide == deriv % diff_nuclide) then
! Find the index in this material for the diff_nuclide.
do j = 1, mat % n_nuclides
@ -3722,7 +3722,7 @@ contains
case (DIFF_TEMPERATURE)
associate (mat => materials(p % material))
if (mat % id == deriv % diff_material) then
if (mat % id() == deriv % diff_material) then
do l=1, mat % n_nuclides
associate (nuc => nuclides(mat % nuclide(l)))
if (mat % nuclide(l) == p % event_nuclide .and. &

View file

@ -76,7 +76,7 @@ contains
allocate(material_ids(size(this % materials)))
do i = 1, size(this % materials)
material_ids(i) = materials(this % materials(i)) % id
material_ids(i) = materials(this % materials(i)) % id()
end do
call write_dataset(filter_group, "bins", material_ids)
end subroutine to_statepoint_material
@ -110,7 +110,7 @@ contains
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
label = "Material " // to_str(materials(this % materials(bin)) % id)
label = "Material " // to_str(materials(this % materials(bin)) % id())
end function text_label_material
!===============================================================================

View file

@ -196,7 +196,7 @@ contains
i_material = p % material
if (i_material /= MATERIAL_VOID) then
do i_domain = 1, size(this % domain_id)
if (materials(i_material) % id == this % domain_id(i_domain)) then
if (materials(i_material) % id() == this % domain_id(i_domain)) then
call check_hit(i_domain, i_material, indices, hits, n_mat)
end if
end do

View file

@ -17,13 +17,14 @@ check_for_node(pugi::xml_node node, const char *name)
}
std::string get_node_value(pugi::xml_node node, const char *name,
bool lowercase=false, bool strip=false);
bool lowercase=false, bool strip=false);
template <typename T>
std::vector<T> get_node_array(pugi::xml_node node, const char* name)
std::vector<T> get_node_array(pugi::xml_node node, const char* name,
bool lowercase=false)
{
// Get value of node attribute/child
std::string s {get_node_value(node, name)};
std::string s {get_node_value(node, name, lowercase)};
// Read values one by one into vector
std::stringstream iss {s};