diff --git a/src/cell.cpp b/src/cell.cpp index 52b9a2e44f..e48a3903a2 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -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 { diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 8399817fa4..7167f24a85 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -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 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::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(); } diff --git a/src/material.cpp b/src/material.cpp index 3bc2452d39..7979fc5f09 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -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) { diff --git a/src/material_header.F90 b/src/material_header.F90 index 3a68a48694..a9de08ec4c 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -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 diff --git a/src/surface.cpp b/src/surface.cpp index 99cfd82b05..b18d7c9195 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -27,7 +27,7 @@ extern "C" const int BC_PERIODIC {3}; int32_t n_surfaces; -Surface** surfaces_c; +std::vector global_surfaces; std::map 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(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(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(); } diff --git a/src/surface.h b/src/surface.h index d43b24ae06..746db5c3b5 100644 --- a/src/surface.h +++ b/src/surface.h @@ -4,6 +4,7 @@ #include #include // For numeric_limits #include +#include #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 global_surfaces; extern std::map surface_map;