Remove cells array on Fortran side

This commit is contained in:
Paul Romano 2019-02-15 06:35:59 -06:00
parent 47c19353e4
commit 964fdfde30
20 changed files with 162 additions and 588 deletions

View file

@ -310,7 +310,6 @@ add_library(libopenmc SHARED
src/dict_header.F90
src/error.F90
src/geometry.F90
src/geometry_header.F90
src/hdf5_interface.F90
src/initialize.F90
src/input_xml.F90

View file

@ -199,6 +199,8 @@ public:
// Non-member functions
//==============================================================================
void read_cells(pugi::xml_node node);
#ifdef DAGMC
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed);
#endif

View file

@ -26,6 +26,7 @@ extern moab::DagMC* DAG;
extern "C" void load_dagmc_geometry();
extern "C" void free_memory_dagmc();
void read_geometry_dagmc();
extern "C" pugi::xml_document* read_uwuw_materials();
bool get_uwuw_materials_xml(std::string& s);

View file

@ -10,6 +10,8 @@
namespace openmc {
void read_geometry_xml();
//==============================================================================
//! Replace Universe, Lattice, and Material IDs with indices.
//==============================================================================
@ -111,7 +113,7 @@ extern "C" int maximum_levels(int32_t univ);
//! Deallocates global vectors and maps for cells, universes, and lattices.
//==============================================================================
extern "C" void free_memory_geometry_c();
extern "C" void free_memory_geometry();
} // namespace openmc
#endif // OPENMC_GEOMETRY_AUX_H

View file

@ -273,5 +273,11 @@ private:
std::array<double, 2> pitch_; //!< Lattice tile width and height
};
//==============================================================================
// Non-member functions
//==============================================================================
void read_lattices(pugi::xml_node node);
} // namespace openmc
#endif // OPENMC_LATTICE_H

View file

@ -395,6 +395,12 @@ public:
void to_hdf5_inner(hid_t group_id) const;
};
//==============================================================================
// Non-member functions
//==============================================================================
void read_surfaces(pugi::xml_node node);
//==============================================================================
// Fortran compatibility functions
//==============================================================================

View file

@ -38,6 +38,7 @@ _dll.openmc_cell_set_temperature.errcheck = _error_handler
_dll.openmc_get_cell_index.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_cell_index.restype = c_int
_dll.openmc_get_cell_index.errcheck = _error_handler
_dll.cells_size.restype = c_int
class Cell(_FortranObjectWithID):
@ -156,10 +157,10 @@ class _CellMapping(Mapping):
def __iter__(self):
for i in range(len(self)):
yield Cell(index=i + 1).id
yield Cell(index=i).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_cells').value
return _dll.cells_size()
def __repr__(self):
return repr(dict(self))

View file

@ -39,7 +39,7 @@ contains
err = E_UNASSIGNED
if (found) then
index = p % coord(p % n_coord) % cell + 1
index = p % coord(p % n_coord) % cell
instance = p % cell_instance
err = 0
else
@ -59,7 +59,6 @@ contains
subroutine free_memory() bind(C)
use bank_header
use geometry_header
use material_header
use nuclide_header
use photon_header
@ -93,6 +92,9 @@ contains
subroutine free_memory_surfaces() bind(C)
end subroutine
subroutine free_memory_geometry() bind(C)
end subroutine
subroutine sab_clear() bind(C)
end subroutine
end interface

View file

@ -23,8 +23,6 @@ namespace openmc {
namespace model {
int32_t n_cells {0};
std::vector<Cell*> cells;
std::unordered_map<int32_t, int32_t> cell_map;
@ -643,18 +641,18 @@ void DAGCell::to_hdf5(hid_t group_id) const { return; }
// Non-method functions
//==============================================================================
extern "C" void
read_cells(pugi::xml_node* node)
void read_cells(pugi::xml_node node)
{
// Count the number of cells.
for (pugi::xml_node cell_node: node->children("cell")) {model::n_cells++;}
if (model::n_cells == 0) {
int n_cells = 0;
for (pugi::xml_node cell_node: node.children("cell")) {n_cells++;}
if (n_cells == 0) {
fatal_error("No cells found in geometry.xml!");
}
// Loop over XML cell elements and populate the array.
model::cells.reserve(model::n_cells);
for (pugi::xml_node cell_node : node->children("cell")) {
model::cells.reserve(n_cells);
for (pugi::xml_node cell_node : node.children("cell")) {
model::cells.push_back(new CSGCell(cell_node));
}
@ -699,9 +697,8 @@ read_cells(pugi::xml_node* node)
extern "C" int
openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
{
if (index >= 1 && index <= model::cells.size()) {
//TODO: off-by-one
Cell& c {*model::cells[index - 1]};
if (index >= 0 && index < model::cells.size()) {
Cell& c {*model::cells[index]};
*type = c.type_;
if (c.type_ == FILL_MATERIAL) {
*indices = c.material_.data();
@ -721,9 +718,8 @@ extern "C" int
openmc_cell_set_fill(int32_t index, int type, int32_t n,
const int32_t* indices)
{
if (index >= 1 && index <= model::cells.size()) {
//TODO: off-by-one
Cell& c {*model::cells[index - 1]};
if (index >= 0 && index < model::cells.size()) {
Cell& c {*model::cells[index]};
if (type == FILL_MATERIAL) {
c.type_ = FILL_MATERIAL;
c.material_.clear();
@ -756,9 +752,9 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
extern "C" int
openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance)
{
if (index >= 1 && index <= model::cells.size()) {
if (index >= 0 && index < model::cells.size()) {
//TODO: off-by-one
Cell& c {*model::cells[index - 1]};
Cell& c {*model::cells[index]};
if (instance) {
if (*instance >= 0 && *instance < c.sqrtkT_.size()) {
@ -781,28 +777,66 @@ openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance)
return 0;
}
//! Return the index in the cells array of a cell with a given ID
extern "C" int
openmc_get_cell_index(int32_t id, int32_t* index)
{
auto it = model::cell_map.find(id);
if (it != model::cell_map.end()) {
*index = it->second;
return 0;
} else {
set_errmsg("No cell exists with ID=" + std::to_string(id) + ".");
return OPENMC_E_INVALID_ID;
}
}
//! Return the ID of a cell
extern "C" int
openmc_cell_get_id(int32_t index, int32_t* id)
{
if (index >= 0 && index < model::cells.size()) {
*id = model::cells[index]->id_;
return 0;
} else {
set_errmsg("Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
//! Set the ID of a cell
extern "C" int
openmc_cell_set_id(int32_t index, int32_t id)
{
if (index >= 0 && index < model::cells.size()) {
model::cells[index]->id_ = id;
model::cell_map[id] = index;
return 0;
} else {
set_errmsg("Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
//! Extend the cells array by n elements
extern "C" int
openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end)
{
if (index_start) *index_start = model::cells.size();
if (index_end) *index_end = model::cells.size() + n - 1;
for (int32_t i = 0; i < n; i++) {
model::cells.push_back(new CSGCell());
}
return 0;
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
Cell* cell_pointer(int32_t cell_ind) {return model::cells[cell_ind];}
int32_t cell_id(Cell* c) {return c->id_;}
void
cell_set_id(Cell* c, int32_t id)
{
c->id_ = id;
// Find the index of this cell and update the cell map.
for (int i = 0; i < model::cells.size(); i++) {
if (model::cells[i] == c) {
model::cell_map[id] = i;
break;
}
}
}
int cells_size() { return model::cells.size(); }
#ifdef DAGMC
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed)
@ -818,33 +852,6 @@ extern "C" {
return cur_cell->dagmc_ptr_->index_by_handle(new_vol);
}
#endif
int32_t cell_universe(Cell* c) {return c->universe_;}
int32_t cell_fill(Cell* c) {return c->fill_;}
int cell_material_size(Cell* c) {return c->material_.size();}
//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;
}
int cell_sqrtkT_size(Cell* c) {return c->sqrtkT_.size();}
double cell_sqrtkT(Cell* c, int i) {return c->sqrtkT_[i];}
void extend_cells_c(int32_t n)
{
model::cells.reserve(model::cells.size() + n);
for (int32_t i = 0; i < n; i++) {
model::cells.push_back(new CSGCell());
}
model::n_cells = model::cells.size();
}
}

View file

@ -4,6 +4,7 @@
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/geometry_aux.h"
#include "openmc/string_utils.h"
#include "openmc/settings.h"
#include "openmc/geometry.h"
@ -133,9 +134,9 @@ void load_dagmc_geometry()
/// Cells (Volumes) \\\
// initialize cell objects
model::n_cells = model::DAG->num_entities(3);
int n_cells = model::DAG->num_entities(3);
moab::EntityHandle graveyard = 0;
for (int i = 0; i < model::n_cells; i++) {
for (int i = 0; i < n_cells; i++) {
moab::EntityHandle vol_handle = model::DAG->entity_by_index(3, i+1);
// set cell ids using global IDs
@ -312,6 +313,20 @@ void load_dagmc_geometry()
return;
}
void read_geometry_dagmc()
{
// Check if dagmc.h5m exists
std::string filename = settings::path_input + "dagmc.h5m"
if (!file_exists(filename)) {
fatal_error("Geometry DAGMC file '" + filename + "' does not exist!");
}
write_message("Reading DAGMC geometry...", 5);
load_dagmc_geometry();
model::root_universe = find_root_universe()
}
void free_memory_dagmc()
{
delete model::DAG;

View file

@ -1,7 +1,6 @@
module geometry
use constants
use geometry_header
use particle_header
use simulation_header
use settings

View file

@ -4,10 +4,14 @@
#include <sstream>
#include <unordered_set>
#include "pugixml.hpp"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/container_util.h"
#include "openmc/dagmc.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/geometry.h"
#include "openmc/lattice.h"
#include "openmc/material.h"
@ -19,6 +23,46 @@
namespace openmc {
void read_geometry_xml()
{
#ifdef DAGMC
if (settings::dagmc) {
read_geometry_dagmc();
return;
}
#endif
// Display output message
write_message("Reading geometry XML file...", 5);
// Check if geometry.xml exists
std::string filename = settings::path_input + "geometry.xml";
if (!file_exists(filename)) {
fatal_error("Geometry XML file '" + filename + "' does not exist!");
}
// Parse settings.xml file
pugi::xml_document doc;
auto result = doc.load_file(filename.c_str());
if (!result) {
fatal_error("Error processing geometry.xml file.");
}
// Get root element
pugi::xml_node root = doc.document_element();
// Read surfaces, cells, lattice
read_surfaces(root);
read_cells(root);
read_lattices(root);
// ==========================================================================
// SETUP UNIVERSES
// Allocate universes, universe cell arrays, and assign base universe
model::root_universe = find_root_universe();
}
//==============================================================================
void
@ -473,12 +517,11 @@ maximum_levels(int32_t univ)
//==============================================================================
void
free_memory_geometry_c()
free_memory_geometry()
{
for (Cell* c : model::cells) {delete c;}
model::cells.clear();
model::cell_map.clear();
model::n_cells = 0;
for (Universe* u : model::universes) {delete u;}
model::universes.clear();

View file

@ -1,279 +0,0 @@
module geometry_header
use, intrinsic :: ISO_C_BINDING
use dict_header, only: DictIntInt
use error
use string, only: to_str
implicit none
interface
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
function cell_id_c(cell_ptr) bind(C, name='cell_id') result(id)
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT32_T) :: id
end function cell_id_c
subroutine cell_set_id_c(cell_ptr, id) bind(C, name='cell_set_id')
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT32_T), intent(in), value :: id
end subroutine cell_set_id_c
function cell_universe_c(cell_ptr) bind(C, name='cell_universe') &
result(universe)
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT32_T) :: universe
end function cell_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_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_sqrtkT_size_c(cell_ptr) bind(C, name='cell_sqrtkT_size') &
result(n)
import C_PTR, C_INT
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT) :: n
end function cell_sqrtkT_size_c
function cell_sqrtkT_c(cell_ptr, i) bind(C, name='cell_sqrtkT') &
result(sqrtkT)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT), intent(in), value :: i
real(C_DOUBLE) :: sqrtkT
end function cell_sqrtkT_c
subroutine extend_cells_c(n) bind(C)
import C_INT32_t
integer(C_INT32_T), intent(in), value :: n
end subroutine extend_cells_c
end interface
!===============================================================================
! CELL defines a closed volume by its bounding surfaces
!===============================================================================
type Cell
type(C_PTR) :: ptr
contains
procedure :: id => cell_id
procedure :: set_id => cell_set_id
procedure :: universe => cell_universe
procedure :: fill => cell_fill
procedure :: material_size => cell_material_size
procedure :: material => cell_material
procedure :: sqrtkT_size => cell_sqrtkT_size
procedure :: sqrtkT => cell_sqrtkT
end type Cell
! array index of the root universe
integer(C_INT), bind(C) :: root_universe
integer(C_INT32_T), bind(C) :: n_cells ! # of cells
integer(C_INT32_T), bind(C) :: n_universes ! # of universes
type(Cell), allocatable, target :: cells(:)
! Dictionaries which map user IDs to indices in the global arrays
type(DictIntInt) :: cell_dict
contains
!===============================================================================
function cell_id(this) result(id)
class(Cell), intent(in) :: this
integer(C_INT32_T) :: id
id = cell_id_c(this % ptr)
end function cell_id
subroutine cell_set_id(this, id)
class(Cell), intent(in) :: this
integer(C_INT32_T), intent(in) :: id
call cell_set_id_c(this % ptr, id)
end subroutine cell_set_id
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
function cell_fill(this) result(fill)
class(Cell), intent(in) :: this
integer(C_INT32_T) :: fill
fill = cell_fill_c(this % ptr)
end function cell_fill
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_sqrtkT_size(this) result(n)
class(Cell), intent(in) :: this
integer :: n
n = cell_sqrtkT_size_c(this % ptr)
end function cell_sqrtkT_size
function cell_sqrtkT(this, i) result(sqrtkT)
class(Cell), intent(in) :: this
integer, intent(in) :: i
real(C_DOUBLE) :: sqrtkT
sqrtkT = cell_sqrtkT_c(this % ptr, i)
end function cell_sqrtkT
!===============================================================================
! FREE_MEMORY_GEOMETRY deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_geometry()
interface
subroutine free_memory_geometry_c() bind(C)
end subroutine free_memory_geometry_c
end interface
call free_memory_geometry_c()
n_cells = 0
n_universes = 0
if (allocated(cells)) deallocate(cells)
call cell_dict % clear()
end subroutine free_memory_geometry
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_extend_cells(n, index_start, index_end) result(err) bind(C)
! Extend the cells array by n elements
integer(C_INT32_T), value, intent(in) :: n
integer(C_INT32_T), optional, intent(out) :: index_start
integer(C_INT32_T), optional, intent(out) :: index_end
integer(C_INT32_T) :: i
integer(C_INT) :: err
type(Cell), allocatable :: temp(:) ! temporary cells array
if (n_cells == 0) then
! Allocate cells array
allocate(cells(n))
else
! Allocate cells array with increased size
allocate(temp(n_cells + n))
! Copy original cells to temporary array
temp(1:n_cells) = cells
! Move allocation from temporary array
call move_alloc(FROM=temp, TO=cells)
end if
! Return indices in cells array
if (present(index_start)) index_start = n_cells + 1
if (present(index_end)) index_end = n_cells + n
! 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(i - 1)
end do
err = 0
end function openmc_extend_cells
function openmc_get_cell_index(id, index) result(err) bind(C)
! Return the index in the cells array of a cell with a given ID
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(cells)) then
if (cell_dict % has(id)) then
index = cell_dict % get(id)
err = 0
else
err = E_INVALID_ID
call set_errmsg("No cell exists with ID=" // trim(to_str(id)) // ".")
end if
else
err = E_ALLOCATE
call set_errmsg("Memory has not been allocated for cells.")
end if
end function openmc_get_cell_index
function openmc_cell_get_id(index, id) result(err) bind(C)
! Return the ID of a cell
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(cells)) then
id = cells(index) % id()
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in cells array is out of bounds.")
end if
end function openmc_cell_get_id
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
integer(C_INT32_T), value, intent(in) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= n_cells) then
call cells(index) % set_id(id)
call cell_dict % set(id, index)
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in cells array is out of bounds.")
end if
end function openmc_cell_set_id
end module geometry_header

View file

@ -31,7 +31,6 @@
// data/functions from Fortran side
extern "C" void read_command_line();
extern "C" void read_geometry_xml();
extern "C" void read_materials_xml();
extern "C" void read_plots_xml();
extern "C" void read_tallies_xml();

View file

@ -3,9 +3,7 @@ module input_xml
use, intrinsic :: ISO_C_BINDING
use constants
use dict_header, only: DictIntInt, DictCharInt, DictEntryCI
use error, only: fatal_error, warning, write_message, openmc_err_msg
use geometry_header
#ifdef DAGMC
use dagmc_header
#endif
@ -30,42 +28,11 @@ module input_xml
save
interface
subroutine count_cell_instances(univ_indx) bind(C)
import C_INT32_T
integer(C_INT32_T), intent(in), value :: univ_indx
end subroutine count_cell_instances
subroutine read_surfaces(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine read_surfaces
subroutine read_cells(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine read_cells
subroutine read_lattices(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine read_lattices
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
end function find_root_universe
function maximum_levels(univ) bind(C) result(n)
import C_INT32_T, C_INT
integer(C_INT32_T), intent(in), value :: univ
integer(C_INT) :: n
end function maximum_levels
subroutine read_plots(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
@ -81,199 +48,6 @@ module input_xml
contains
#ifdef DAGMC
!===============================================================================
! READ_GEOMETRY_DAGMC reads data from a DAGMC .h5m file, checking
! for material properties and surface boundary conditions
! some universe information is spoofed for now
!===============================================================================
subroutine read_geometry_dagmc()
integer :: i
integer :: univ_id
integer :: n_cells_in_univ
logical :: file_exists
character(MAX_LINE_LEN) :: filename
type(Cell), pointer :: c
type(VectorInt) :: univ_ids ! List of all universe IDs
type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
! universe contains
! Check if dagmc.h5m exists
filename = trim(path_input) // "dagmc.h5m"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Geometry DAGMC file '" // trim(filename) // "' does not &
&exist!")
end if
call write_message("Reading DAGMC geometry...", 5)
call load_dagmc_geometry()
call allocate_cells()
! setup universe data structs
do i = 1, n_cells
c => cells(i)
! additional metadata spoofing
univ_id = c % universe()
if (.not. cells_in_univ_dict % has(univ_id)) then
n_universes = n_universes + 1
n_cells_in_univ = 1
call univ_ids % push_back(univ_id)
else
n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id)
end if
call cells_in_univ_dict % set(univ_id, n_cells_in_univ)
end do
root_universe = find_root_universe()
end subroutine read_geometry_dagmc
#endif
!===============================================================================
! READ_GEOMETRY_XML reads data from a geometry.xml file and parses it, checking
! for errors and placing properly-formatted data in the right data structures
!===============================================================================
subroutine read_geometry_xml() bind(C)
integer :: i, n
integer :: univ_id
integer :: n_cells_in_univ
real(8) :: phi, theta, psi
logical :: file_exists
character(MAX_LINE_LEN) :: filename
type(Cell), pointer :: c
type(XMLDocument) :: doc
type(XMLNode) :: root
type(XMLNode) :: node_cell
type(XMLNode), allocatable :: node_cell_list(:)
type(VectorInt) :: univ_ids ! List of all universe IDs
type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
! universe contains
#ifdef DAGMC
if (dagmc) then
call read_geometry_dagmc()
return
end if
#endif
! Display output message
call write_message("Reading geometry XML file...", 5)
! Check if geometry.xml exists
filename = trim(path_input) // "geometry.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Geometry XML file '" // trim(filename) // "' does not &
&exist!")
end if
! Parse geometry.xml file
call doc % load_file(filename)
root = doc % document_element()
! ==========================================================================
! READ SURFACES FROM GEOMETRY.XML
call read_surfaces(root % ptr)
! ==========================================================================
! READ CELLS FROM GEOMETRY.XML
call read_cells(root % ptr)
! Get pointer to list of XML <cell>
call get_node_list(root, "cell", node_cell_list)
! Get number of <cell> tags
n_cells = size(node_cell_list)
! Check for no cells
if (n_cells == 0) then
call fatal_error("No cells found in geometry.xml!")
end if
! Allocate cells array
allocate(cells(n_cells))
n_universes = 0
do i = 1, n_cells
c => cells(i)
c % ptr = cell_pointer(i - 1)
! Get pointer to i-th cell node
node_cell = node_cell_list(i)
! Check to make sure 'id' hasn't been used
if (cell_dict % has(c % id())) then
call fatal_error("Two or more cells use the same unique ID: " &
// to_str(c % id()))
end if
! Add cell to dictionary
call cell_dict % set(c % id(), i)
! For cells, we also need to check if there's a new universe --
! also for every cell add 1 to the count of cells for the
! specified universe
univ_id = c % universe()
if (.not. cells_in_univ_dict % has(univ_id)) then
n_universes = n_universes + 1
n_cells_in_univ = 1
call univ_ids % push_back(univ_id)
else
n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id)
end if
call cells_in_univ_dict % set(univ_id, n_cells_in_univ)
end do
! ==========================================================================
! READ LATTICES FROM GEOMETRY.XML
call read_lattices(root % ptr)
! ==========================================================================
! SETUP UNIVERSES
! Allocate universes, universe cell arrays, and assign base universe
root_universe = find_root_universe()
! Clear dictionary
call cells_in_univ_dict%clear()
! Close geometry XML file
call doc % clear()
end subroutine read_geometry_xml
subroutine allocate_cells()
integer :: i
type(Cell), pointer :: c
! Allocate cells array
allocate(cells(n_cells))
do i = 1, n_cells
c => cells(i)
c % ptr = cell_pointer(i - 1)
! Check to make sure 'id' hasn't been used
if (cell_dict % has(c % id())) then
call fatal_error("Two or more cells use the same unique ID: " &
// to_str(c % id()))
end if
! Add cell to dictionary
call cell_dict % set(c % id(), i)
end do
end subroutine allocate_cells
subroutine read_materials_xml() bind(C)
logical :: file_exists ! does materials.xml exist?
character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml

View file

@ -864,13 +864,12 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const
// Non-method functions
//==============================================================================
extern "C" void
read_lattices(pugi::xml_node *node)
void read_lattices(pugi::xml_node node)
{
for (pugi::xml_node lat_node : node->children("lattice")) {
for (pugi::xml_node lat_node : node.children("lattice")) {
model::lattices.push_back(new RectLattice(lat_node));
}
for (pugi::xml_node lat_node : node->children("hex_lattice")) {
for (pugi::xml_node lat_node : node.children("hex_lattice")) {
model::lattices.push_back(new HexLattice(lat_node));
}

View file

@ -172,7 +172,7 @@ Bank SourceDistribution::sample() const
if (space_box) {
if (space_box->only_fissionable()) {
// Determine material
auto c = model::cells[cell_index - 1];
auto c = model::cells[cell_index];
auto mat_index = c->material_.size() == 1
? c->material_[0] : c->material_[instance];

View file

@ -1061,12 +1061,11 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
//==============================================================================
extern "C" void
read_surfaces(pugi::xml_node* node)
void read_surfaces(pugi::xml_node node)
{
// Count the number of surfaces.
int n_surfaces = 0;
for (pugi::xml_node surf_node : node->children("surface")) {n_surfaces++;}
for (pugi::xml_node surf_node : node.children("surface")) {n_surfaces++;}
if (n_surfaces == 0) {
fatal_error("No surfaces found in geometry.xml!");
}
@ -1076,7 +1075,7 @@ read_surfaces(pugi::xml_node* node)
{
pugi::xml_node surf_node;
int i_surf;
for (surf_node = node->child("surface"), i_surf = 0; surf_node;
for (surf_node = node.child("surface"), i_surf = 0; surf_node;
surf_node = surf_node.next_sibling("surface"), i_surf++) {
std::string surf_type = get_node_value(surf_node, "type", true, true);

View file

@ -4,7 +4,6 @@ module tally
use bank_header
use constants
use geometry_header
use material_header
use message_passing
use mgxs_interface

View file

@ -72,7 +72,7 @@ def test_cell(capi_init):
cell = openmc.capi.cells[1]
assert isinstance(cell.fill, openmc.capi.Material)
cell.fill = openmc.capi.materials[1]
assert str(cell) == 'Cell[1]'
assert str(cell) == 'Cell[0]'
def test_new_cell(capi_init):