This commit is contained in:
Sterling Harper 2018-08-15 21:41:09 -04:00
parent 1f11eabc71
commit 20eeb75c90
6 changed files with 80 additions and 71 deletions

View file

@ -204,9 +204,8 @@ Cell::Cell(pugi::xml_node cell_node)
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");
name = get_node_value(cell_node, "name", false);
}
if (check_for_node(cell_node, "universe")) {
@ -238,8 +237,8 @@ Cell::Cell(pugi::xml_node cell_node)
}
} else {
material.push_back(C_NONE);
material.shrink_to_fit();
}
material.shrink_to_fit();
} else {
material.push_back(C_NONE);
material.shrink_to_fit();
@ -318,7 +317,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) {
@ -359,7 +358,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);
@ -382,7 +382,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;}
}
}
@ -424,7 +424,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));
}
}
@ -454,13 +454,12 @@ 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));
}
global_cells.shrink_to_fit();
// Populate the Universe vector and map.
for (int i = 0; i < global_cells.size(); i++) {
@ -475,6 +474,7 @@ read_cells(pugi::xml_node* node)
global_universes[it->second]->cells.push_back(i);
}
}
global_universes.shrink_to_fit();
}
//==============================================================================
@ -524,6 +524,7 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
return OPENMC_E_OUT_OF_BOUNDS;
}
}
c.material.shrink_to_fit();
} else if (type == FILL_UNIVERSE) {
c.type = FILL_UNIVERSE;
} else {

View file

@ -19,7 +19,7 @@ void
adjust_indices()
{
// Adjust material/fill idices.
for (Cell *c : global_cells) {
for (Cell* c : global_cells) {
if (c->material[0] == C_NONE) {
int32_t id = c->fill;
auto search_univ = universe_map.find(id);
@ -56,7 +56,7 @@ adjust_indices()
}
// 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.
@ -70,7 +70,7 @@ adjust_indices()
}
// Change all lattice universe values from IDs to indices.
for (Lattice *l : lattices_c) {
for (Lattice* l : lattices_c) {
l->adjust_indices();
}
}
@ -82,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);
}
@ -123,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);
}
}
@ -140,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) {
@ -149,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);
}
@ -169,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);
@ -192,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;
@ -203,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);
}
}
@ -214,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;
@ -224,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();
}
@ -236,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) {
@ -244,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];
}
@ -256,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) {
@ -267,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;
@ -289,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;
@ -299,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());
@ -316,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));
@ -338,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

@ -43,6 +43,7 @@ read_materials(pugi::xml_node* node)
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++) {
@ -67,7 +68,12 @@ extern "C" {
int32_t material_id(Material* mat) {return mat->id;}
void material_set_id(Material* mat, int32_t id) {mat->id = 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)
{

View file

@ -42,10 +42,12 @@ module material_header
integer(C_INT32_T) :: id
end function material_id_c
subroutine material_set_id_c(mat_ptr, id) bind(C, name='material_set_id')
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)
@ -122,10 +124,11 @@ contains
id = material_id_c(this % ptr)
end function material_id
subroutine material_set_id(this, id)
subroutine material_set_id(this, id, index)
class(Material), intent(in) :: this
integer(C_INT32_T), intent(in) :: id
call material_set_id_c(this % ptr, 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)
@ -674,7 +677,7 @@ contains
integer(C_INT) :: err
if (index >= 1 .and. index <= n_materials) then
call materials(index) % set_id(id)
call materials(index) % set_id(id, index)
call material_dict % set(id, index)
err = 0
else

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;
@ -147,7 +147,7 @@ Surface::Surface(pugi::xml_node surf_node)
}
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")) {
@ -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;
@ -1085,10 +1083,11 @@ read_surfaces(pugi::xml_node* node)
}
}
}
global_surfaces.shrink_to_fit();
// 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,9 +1103,9 @@ 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];
Surface* surf_base = global_surfaces[i_surf];
PeriodicSurface* surf = dynamic_cast<PeriodicSurface*>(surf_base);
// Make sure this surface inherits from PeriodicSurface.
@ -1149,9 +1148,9 @@ 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];
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
@ -1198,7 +1197,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,7 +1212,7 @@ 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;}
@ -1264,9 +1263,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;