mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Move more Surface implementation to C++
This commit is contained in:
Sterling Harper
parent
afe2f619c8
commit
164e6c0ef4
8 changed files with 217 additions and 248 deletions
|
|
@ -54,16 +54,24 @@ module geometry
|
|||
real(C_DOUBLE), intent(out) :: uvw(3);
|
||||
end subroutine surface_normal_c
|
||||
|
||||
function surface_periodic_c(surf_ind1, surf_ind2, xyz, uvw) &
|
||||
function surface_periodic_c(surf_ind1, xyz, uvw) &
|
||||
bind(C, name="surface_periodic") result(rotational)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
integer(C_INT), intent(in), value :: surf_ind1;
|
||||
integer(C_INT), intent(in), value :: surf_ind2;
|
||||
real(C_DOUBLE), intent(inout) :: xyz(3);
|
||||
real(C_DOUBLE), intent(inout) :: uvw(3);
|
||||
logical(C_BOOL) :: rotational
|
||||
end function surface_periodic_c
|
||||
|
||||
function surface_i_periodic_c(surf_ind) bind(C, name="surface_i_periodic") &
|
||||
result(i_periodic)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
integer(C_INT), intent(in), value :: surf_ind
|
||||
integer(C_INT) :: i_periodic
|
||||
end function
|
||||
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
|
@ -857,15 +865,15 @@ contains
|
|||
! Copy positive neighbors to Surface instance
|
||||
n = neighbor_pos(i)%size()
|
||||
if (n > 0) then
|
||||
allocate(surfaces(i)%obj%neighbor_pos(n))
|
||||
surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:n)
|
||||
allocate(surfaces(i)%neighbor_pos(n))
|
||||
surfaces(i)%neighbor_pos(:) = neighbor_pos(i)%data(1:n)
|
||||
end if
|
||||
|
||||
! Copy negative neighbors to Surface instance
|
||||
n = neighbor_neg(i)%size()
|
||||
if (n > 0) then
|
||||
allocate(surfaces(i)%obj%neighbor_neg(n))
|
||||
surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:n)
|
||||
allocate(surfaces(i)%neighbor_neg(n))
|
||||
surfaces(i)%neighbor_neg(:) = neighbor_neg(i)%data(1:n)
|
||||
end if
|
||||
end do
|
||||
|
||||
|
|
|
|||
|
|
@ -919,12 +919,8 @@ contains
|
|||
integer :: id
|
||||
integer :: univ_id
|
||||
integer :: n_cells_in_univ
|
||||
integer :: coeffs_reqd
|
||||
integer :: i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax
|
||||
real(8) :: xmin, xmax, ymin, ymax, zmin, zmax
|
||||
integer, allocatable :: temp_int_array(:)
|
||||
real(8) :: phi, theta, psi
|
||||
real(8), allocatable :: coeffs(:)
|
||||
logical :: file_exists
|
||||
logical :: boundary_exists
|
||||
character(MAX_LINE_LEN) :: filename
|
||||
|
|
@ -984,13 +980,6 @@ contains
|
|||
call fatal_error("No surfaces found in geometry.xml!")
|
||||
end if
|
||||
|
||||
xmin = INFINITY
|
||||
xmax = -INFINITY
|
||||
ymin = INFINITY
|
||||
ymax = -INFINITY
|
||||
zmin = INFINITY
|
||||
zmax = -INFINITY
|
||||
|
||||
! Allocate cells array
|
||||
allocate(surfaces(n_surfaces))
|
||||
|
||||
|
|
@ -998,52 +987,7 @@ contains
|
|||
! Get pointer to i-th surface node
|
||||
node_surf = node_surf_list(i)
|
||||
|
||||
! Copy and interpret surface type
|
||||
word = ''
|
||||
if (check_for_node(node_surf, "type")) &
|
||||
call get_node_value(node_surf, "type", word)
|
||||
select case(to_lower(word))
|
||||
case ('x-plane')
|
||||
coeffs_reqd = 1
|
||||
allocate(SurfaceXPlane :: surfaces(i)%obj)
|
||||
case ('y-plane')
|
||||
coeffs_reqd = 1
|
||||
allocate(SurfaceYPlane :: surfaces(i)%obj)
|
||||
case ('z-plane')
|
||||
coeffs_reqd = 1
|
||||
allocate(SurfaceZPlane :: surfaces(i)%obj)
|
||||
case ('plane')
|
||||
coeffs_reqd = 4
|
||||
allocate(SurfacePlane :: surfaces(i)%obj)
|
||||
case ('x-cylinder')
|
||||
coeffs_reqd = 3
|
||||
allocate(SurfaceXCylinder :: surfaces(i)%obj)
|
||||
case ('y-cylinder')
|
||||
coeffs_reqd = 3
|
||||
allocate(SurfaceYCylinder :: surfaces(i)%obj)
|
||||
case ('z-cylinder')
|
||||
coeffs_reqd = 3
|
||||
allocate(SurfaceZCylinder :: surfaces(i)%obj)
|
||||
case ('sphere')
|
||||
coeffs_reqd = 4
|
||||
allocate(SurfaceSphere :: surfaces(i)%obj)
|
||||
case ('x-cone')
|
||||
coeffs_reqd = 4
|
||||
allocate(SurfaceXCone :: surfaces(i)%obj)
|
||||
case ('y-cone')
|
||||
coeffs_reqd = 4
|
||||
allocate(SurfaceYCone :: surfaces(i)%obj)
|
||||
case ('z-cone')
|
||||
coeffs_reqd = 4
|
||||
allocate(SurfaceZCone :: surfaces(i)%obj)
|
||||
case ('quadric')
|
||||
coeffs_reqd = 10
|
||||
allocate(SurfaceQuadric :: surfaces(i)%obj)
|
||||
case default
|
||||
call fatal_error("Invalid surface type: " // trim(word))
|
||||
end select
|
||||
|
||||
s => surfaces(i)%obj
|
||||
s => surfaces(i)
|
||||
|
||||
! Copy data into cells
|
||||
if (check_for_node(node_surf, "id")) then
|
||||
|
|
@ -1058,44 +1002,10 @@ contains
|
|||
// to_str(s%id))
|
||||
end if
|
||||
|
||||
! Copy surface name
|
||||
if (check_for_node(node_surf, "name")) then
|
||||
call get_node_value(node_surf, "name", s%name)
|
||||
end if
|
||||
|
||||
! Check to make sure that the proper number of coefficients
|
||||
! have been specified for the given type of surface. Then copy
|
||||
! surface coordinates.
|
||||
|
||||
n = node_word_count(node_surf, "coeffs")
|
||||
|
||||
allocate(coeffs(n))
|
||||
call get_node_array(node_surf, "coeffs", coeffs)
|
||||
|
||||
select type(s)
|
||||
type is (SurfaceXPlane)
|
||||
! Determine outer surfaces
|
||||
xmin = min(xmin, coeffs(1))
|
||||
xmax = max(xmax, coeffs(1))
|
||||
if (xmin == coeffs(1)) i_xmin = i
|
||||
if (xmax == coeffs(1)) i_xmax = i
|
||||
type is (SurfaceYPlane)
|
||||
! Determine outer surfaces
|
||||
ymin = min(ymin, coeffs(1))
|
||||
ymax = max(ymax, coeffs(1))
|
||||
if (ymin == coeffs(1)) i_ymin = i
|
||||
if (ymax == coeffs(1)) i_ymax = i
|
||||
type is (SurfaceZPlane)
|
||||
! Determine outer surfaces
|
||||
zmin = min(zmin, coeffs(1))
|
||||
zmax = max(zmax, coeffs(1))
|
||||
if (zmin == coeffs(1)) i_zmin = i
|
||||
if (zmax == coeffs(1)) i_zmax = i
|
||||
end select
|
||||
|
||||
! No longer need coefficients
|
||||
deallocate(coeffs)
|
||||
|
||||
! Boundary conditions
|
||||
word = ''
|
||||
if (check_for_node(node_surf, "boundary")) &
|
||||
|
|
@ -1112,12 +1022,6 @@ contains
|
|||
case ('periodic')
|
||||
s%bc = BC_PERIODIC
|
||||
boundary_exists = .true.
|
||||
|
||||
! Check for specification of periodic surface
|
||||
if (check_for_node(node_surf, "periodic_surface_id")) then
|
||||
call get_node_value(node_surf, "periodic_surface_id", &
|
||||
s % i_periodic)
|
||||
end if
|
||||
case default
|
||||
call fatal_error("Unknown boundary condition '" // trim(word) // &
|
||||
&"' specified on surface " // trim(to_str(s%id)))
|
||||
|
|
@ -1134,76 +1038,6 @@ contains
|
|||
end if
|
||||
end if
|
||||
|
||||
! Determine opposite side for periodic boundaries
|
||||
do i = 1, size(surfaces)
|
||||
if (surfaces(i) % obj % bc == BC_PERIODIC) then
|
||||
select type (surf => surfaces(i) % obj)
|
||||
type is (SurfaceXPlane)
|
||||
if (surf % i_periodic == NONE) then
|
||||
if (i == i_xmin) then
|
||||
surf % i_periodic = i_xmax
|
||||
elseif (i == i_xmax) then
|
||||
surf % i_periodic = i_xmin
|
||||
else
|
||||
call fatal_error("Periodic boundary condition applied to &
|
||||
&interior surface.")
|
||||
end if
|
||||
else
|
||||
surf % i_periodic = surface_dict % get(surf % i_periodic)
|
||||
end if
|
||||
|
||||
type is (SurfaceYPlane)
|
||||
if (surf % i_periodic == NONE) then
|
||||
if (i == i_ymin) then
|
||||
surf % i_periodic = i_ymax
|
||||
elseif (i == i_ymax) then
|
||||
surf % i_periodic = i_ymin
|
||||
else
|
||||
call fatal_error("Periodic boundary condition applied to &
|
||||
&interior surface.")
|
||||
end if
|
||||
else
|
||||
surf % i_periodic = surface_dict % get(surf % i_periodic)
|
||||
end if
|
||||
|
||||
type is (SurfaceZPlane)
|
||||
if (surf % i_periodic == NONE) then
|
||||
if (i == i_zmin) then
|
||||
surf % i_periodic = i_zmax
|
||||
elseif (i == i_zmax) then
|
||||
surf % i_periodic = i_zmin
|
||||
else
|
||||
call fatal_error("Periodic boundary condition applied to &
|
||||
&interior surface.")
|
||||
end if
|
||||
else
|
||||
surf % i_periodic = surface_dict % get(surf % i_periodic)
|
||||
end if
|
||||
|
||||
type is (SurfacePlane)
|
||||
if (surf % i_periodic == NONE) then
|
||||
call fatal_error("No matching periodic surface specified for &
|
||||
&periodic boundary condition on surface " // &
|
||||
trim(to_str(surf % id)) // ".")
|
||||
else
|
||||
surf % i_periodic = surface_dict % get(surf % i_periodic)
|
||||
end if
|
||||
|
||||
class default
|
||||
call fatal_error("Periodic boundary condition applied to &
|
||||
&non-planar surface.")
|
||||
end select
|
||||
|
||||
! Make sure opposite surface is also periodic
|
||||
associate (surf => surfaces(i) % obj)
|
||||
if (surfaces(surf % i_periodic) % obj % bc /= BC_PERIODIC) then
|
||||
call fatal_error("Could not find matching surface for periodic &
|
||||
&boundary on surface " // trim(to_str(surf % id)) // ".")
|
||||
end if
|
||||
end associate
|
||||
end if
|
||||
end do
|
||||
|
||||
! ==========================================================================
|
||||
! READ CELLS FROM GEOMETRY.XML
|
||||
|
||||
|
|
|
|||
|
|
@ -259,7 +259,7 @@ contains
|
|||
|
||||
! Print surface
|
||||
if (p % surface /= NONE) then
|
||||
write(ou,*) ' Surface = ' // to_str(sign(surfaces(i)%obj%id, p % surface))
|
||||
write(ou,*) ' Surface = ' // to_str(sign(surfaces(i)%id, p % surface))
|
||||
end if
|
||||
|
||||
! Display weight, energy, grid index, and interpolation factor
|
||||
|
|
|
|||
|
|
@ -133,12 +133,11 @@ contains
|
|||
real(8), allocatable :: cell_temperatures(:)
|
||||
integer(HID_T) :: geom_group
|
||||
integer(HID_T) :: cells_group, cell_group
|
||||
integer(HID_T) :: surfaces_group, surface_group
|
||||
integer(HID_T) :: surfaces_group
|
||||
integer(HID_T) :: universes_group, univ_group
|
||||
integer(HID_T) :: lattices_group, lattice_group
|
||||
character(:), allocatable :: region_spec
|
||||
type(Cell), pointer :: c
|
||||
class(Surface), pointer :: s
|
||||
class(Lattice), pointer :: lat
|
||||
|
||||
! Use H5LT interface to write number of geometry objects
|
||||
|
|
@ -233,7 +232,7 @@ contains
|
|||
region_spec = trim(region_spec) // " |"
|
||||
case default
|
||||
region_spec = trim(region_spec) // " " // to_str(&
|
||||
sign(surfaces(abs(k))%obj%id, k))
|
||||
sign(surfaces(abs(k))%id, k))
|
||||
end select
|
||||
end do
|
||||
call write_dataset(cell_group, "region", adjustl(region_spec))
|
||||
|
|
|
|||
|
|
@ -1,7 +1,8 @@
|
|||
#include <algorithm> // for std::transform
|
||||
#include <array> // For std::array
|
||||
#include <array>
|
||||
#include <cstring> // For strcmp
|
||||
#include <limits> // For numeric_limits
|
||||
#include <map>
|
||||
#include <math.h> // For fabs
|
||||
|
||||
#include "pugixml/pugixml.hpp"
|
||||
|
|
@ -45,6 +46,7 @@ get_node_value(const pugi::xml_node &node, const char *name)
|
|||
|
||||
std::string value(value_char);
|
||||
std::transform(value.begin(), value.end(), value.begin(), ::tolower);
|
||||
//TODO: trim whitespace
|
||||
return value;
|
||||
}
|
||||
|
||||
|
|
@ -57,6 +59,7 @@ get_node_value(const pugi::xml_node &node, const char *name)
|
|||
extern "C" double FP_COINCIDENT;
|
||||
//extern "C" double INFTY;
|
||||
const double INFTY{std::numeric_limits<double>::max()};
|
||||
const int C_NONE{-1};
|
||||
|
||||
extern "C" const int BC_TRANSMIT{0};
|
||||
extern "C" const int BC_VACUUM{1};
|
||||
|
|
@ -70,6 +73,8 @@ extern "C" const int BC_PERIODIC{3};
|
|||
class Surface;
|
||||
Surface **surfaces_c;
|
||||
|
||||
std::map<int, int> surface_dict;
|
||||
|
||||
//==============================================================================
|
||||
// Helper functions for reading the "coeffs" node of an XML surface element
|
||||
//==============================================================================
|
||||
|
|
@ -182,6 +187,19 @@ void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1, double &c2,
|
|||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
//==============================================================================
|
||||
|
||||
struct BoundingBox
|
||||
{
|
||||
double xmin;
|
||||
double xmax;
|
||||
double ymin;
|
||||
double ymax;
|
||||
double zmin;
|
||||
double zmax;
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
//! A geometry primitive used to define regions of 3D space.
|
||||
//==============================================================================
|
||||
|
|
@ -193,8 +211,6 @@ public:
|
|||
//int neighbor_pos[], //!< List of cells on positive side
|
||||
// neighbor_neg[]; //!< List of cells on negative side
|
||||
int bc; //!< Boundary condition
|
||||
//TODO: switch that zero to a NONE constant.
|
||||
//int i_periodic = 0; //!< Index of corresponding periodic surface
|
||||
std::string name{""}; //!< User-defined name
|
||||
|
||||
Surface(pugi::xml_node surf_node);
|
||||
|
|
@ -249,7 +265,6 @@ Surface::Surface(pugi::xml_node surf_node)
|
|||
} else {
|
||||
fatal_error("Must specify id of surface in geometry XML file.");
|
||||
}
|
||||
//TODO: check for duplicate IDs
|
||||
|
||||
if (check_for_node(surf_node, "name")) {
|
||||
name = get_node_value(surf_node, "name");
|
||||
|
|
@ -362,7 +377,9 @@ Surface::to_hdf5(hid_t group_id) const
|
|||
class PeriodicSurface : public Surface
|
||||
{
|
||||
public:
|
||||
PeriodicSurface(pugi::xml_node surf_node) : Surface(surf_node) {}
|
||||
int i_periodic{C_NONE}; //!< Index of corresponding periodic surface
|
||||
|
||||
PeriodicSurface(pugi::xml_node surf_node);
|
||||
|
||||
//! Translate a particle onto this surface from a periodic partner surface.
|
||||
//! @param other A pointer to the partner surface in this periodic BC.
|
||||
|
|
@ -374,8 +391,18 @@ public:
|
|||
//! boundary condition.
|
||||
virtual bool periodic_translate(PeriodicSurface *other, double xyz[3],
|
||||
double uvw[3]) const = 0;
|
||||
|
||||
virtual struct BoundingBox bounding_box() const = 0;
|
||||
};
|
||||
|
||||
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"));
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Generic functions for x-, y-, and z-, planes.
|
||||
//==============================================================================
|
||||
|
|
@ -428,6 +455,7 @@ public:
|
|||
void to_hdf5_inner(hid_t group_id) const;
|
||||
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
|
||||
const;
|
||||
struct BoundingBox bounding_box() const;
|
||||
};
|
||||
|
||||
SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node)
|
||||
|
|
@ -485,6 +513,13 @@ bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
|
|||
}
|
||||
}
|
||||
|
||||
struct BoundingBox
|
||||
SurfaceXPlane::bounding_box() const
|
||||
{
|
||||
struct BoundingBox out{x0, x0, -INFTY, INFTY, -INFTY, INFTY};
|
||||
return out;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// SurfaceYPlane
|
||||
//! A plane perpendicular to the y-axis.
|
||||
|
|
@ -504,6 +539,7 @@ public:
|
|||
void to_hdf5_inner(hid_t group_id) const;
|
||||
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
|
||||
const;
|
||||
struct BoundingBox bounding_box() const;
|
||||
};
|
||||
|
||||
SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node)
|
||||
|
|
@ -561,6 +597,13 @@ bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
|
|||
}
|
||||
}
|
||||
|
||||
struct BoundingBox
|
||||
SurfaceYPlane::bounding_box() const
|
||||
{
|
||||
struct BoundingBox out{-INFTY, INFTY, y0, y0, -INFTY, INFTY};
|
||||
return out;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// SurfaceZPlane
|
||||
//! A plane perpendicular to the z-axis.
|
||||
|
|
@ -580,6 +623,7 @@ public:
|
|||
void to_hdf5_inner(hid_t group_id) const;
|
||||
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
|
||||
const;
|
||||
struct BoundingBox bounding_box() const;
|
||||
};
|
||||
|
||||
SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node)
|
||||
|
|
@ -619,6 +663,13 @@ bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
|
|||
return false;
|
||||
}
|
||||
|
||||
struct BoundingBox
|
||||
SurfaceZPlane::bounding_box() const
|
||||
{
|
||||
struct BoundingBox out{-INFTY, INFTY, -INFTY, INFTY, z0, z0};
|
||||
return out;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// SurfacePlane
|
||||
//! A general plane.
|
||||
|
|
@ -638,6 +689,7 @@ public:
|
|||
void to_hdf5_inner(hid_t group_id) const;
|
||||
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
|
||||
const;
|
||||
struct BoundingBox bounding_box() const;
|
||||
};
|
||||
|
||||
SurfacePlane::SurfacePlane(pugi::xml_node surf_node)
|
||||
|
|
@ -698,6 +750,13 @@ bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3],
|
|||
return false;
|
||||
}
|
||||
|
||||
struct BoundingBox
|
||||
SurfacePlane::bounding_box() const
|
||||
{
|
||||
struct BoundingBox out{-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
|
||||
return out;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Generic functions for x-, y-, and z-, cylinders
|
||||
//==============================================================================
|
||||
|
|
@ -1359,10 +1418,10 @@ extern "C" void
|
|||
read_surfaces(pugi::xml_node *node)
|
||||
{
|
||||
// Count the number of surfaces.
|
||||
int n_surfaces = 0;
|
||||
int n_surfaces{0};
|
||||
for (pugi::xml_node surf_node = node->child("surface"); surf_node;
|
||||
surf_node = surf_node.next_sibling("surface")) {
|
||||
n_surfaces += 1;
|
||||
n_surfaces++;
|
||||
}
|
||||
|
||||
// Allocate the array of Surface pointers.
|
||||
|
|
@ -1413,12 +1472,125 @@ read_surfaces(pugi::xml_node *node)
|
|||
surfaces_c[i_surf] = new SurfaceQuadric(surf_node);
|
||||
|
||||
} else {
|
||||
std::cout << "Call error or handle uppercase here!" << std::endl;
|
||||
std::cout << "Call error here!" << std::endl;
|
||||
std::cout << surf_type << std::endl;
|
||||
//TODO: call fatal_error
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Fill the surface dictionary.
|
||||
for (int i_surf = 0; i_surf < n_surfaces; i_surf++) {
|
||||
int id = surfaces_c[i_surf]->id;
|
||||
auto in_dict = surface_dict.find(id);
|
||||
if (in_dict == surface_dict.end()) {
|
||||
surface_dict[id] = i_surf;
|
||||
} else {
|
||||
std::string err_msg{"Two or more surfaces use the same unique ID: "};
|
||||
err_msg += std::to_string(id);
|
||||
fatal_error(err_msg);
|
||||
}
|
||||
}
|
||||
|
||||
// Find the global bounding box (of periodic BC surfaces).
|
||||
double xmin{INFTY}, xmax{-INFTY}, ymin{INFTY}, ymax{-INFTY},
|
||||
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) {
|
||||
// Downcast to the PeriodicSurface type.
|
||||
Surface *surf_base = surfaces_c[i_surf];
|
||||
PeriodicSurface *surf = dynamic_cast<PeriodicSurface *>(surf_base);
|
||||
//TODO: check for null pointer
|
||||
|
||||
// See if this surface makes part of the global bounding box.
|
||||
struct BoundingBox bb = surf->bounding_box();
|
||||
if (bb.xmin > -INFTY and bb.xmin < xmin) {
|
||||
xmin = bb.xmin;
|
||||
i_xmin = i_surf;
|
||||
}
|
||||
if (bb.xmax < INFTY and bb.xmax > xmax) {
|
||||
xmax = bb.xmax;
|
||||
i_xmax = i_surf;
|
||||
}
|
||||
if (bb.ymin > -INFTY and bb.ymin < ymin) {
|
||||
ymin = bb.ymin;
|
||||
i_ymin = i_surf;
|
||||
}
|
||||
if (bb.ymax < INFTY and bb.ymax > ymax) {
|
||||
ymax = bb.ymax;
|
||||
i_ymax = i_surf;
|
||||
}
|
||||
if (bb.zmin > -INFTY and bb.zmin < zmin) {
|
||||
zmin = bb.zmin;
|
||||
i_zmin = i_surf;
|
||||
}
|
||||
if (bb.zmax < INFTY and bb.zmax > zmax) {
|
||||
zmax = bb.zmax;
|
||||
i_zmax = i_surf;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 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) {
|
||||
// Downcast to the PeriodicSurface type.
|
||||
Surface *surf_base = surfaces_c[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);
|
||||
|
||||
if (surf_p == NULL) {
|
||||
// This is not a SurfacePlane.
|
||||
if (surf->i_periodic == C_NONE) {
|
||||
// The user did not specify the matching periodic surface. See if we
|
||||
// can find the partnered surface from the bounding box information.
|
||||
if (i_surf == i_xmin) {
|
||||
surf->i_periodic = i_xmax;
|
||||
} else if (i_surf == i_xmax) {
|
||||
surf->i_periodic = i_xmin;
|
||||
} else if (i_surf == i_ymin) {
|
||||
surf->i_periodic = i_ymax;
|
||||
} else if (i_surf == i_ymax) {
|
||||
surf->i_periodic = i_ymin;
|
||||
} else if (i_surf == i_zmin) {
|
||||
surf->i_periodic = i_zmax;
|
||||
} else if (i_surf == i_zmax) {
|
||||
surf->i_periodic = i_zmin;
|
||||
} else {
|
||||
fatal_error("Periodic boundary condition applied to interior "
|
||||
"surface");
|
||||
}
|
||||
} else {
|
||||
// Convert the surface id to an index.
|
||||
surf->i_periodic = surface_dict[surf->i_periodic];
|
||||
}
|
||||
} else {
|
||||
// This is a SurfacePlane. We won't try to find it's partner if the
|
||||
// user didn't specify one.
|
||||
if (surf->i_periodic == C_NONE) {
|
||||
std::string err_msg{"No matching periodic surface specified for "
|
||||
"periodic boundary condition on surface "};
|
||||
err_msg += std::to_string(surf->id);
|
||||
fatal_error(err_msg);
|
||||
} else {
|
||||
// Convert the surface id to an index.
|
||||
surf->i_periodic = surface_dict[surf->i_periodic];
|
||||
}
|
||||
}
|
||||
|
||||
// Make sure the opposite surface is also periodic.
|
||||
if (surfaces_c[surf->i_periodic]->bc != BC_PERIODIC) {
|
||||
std::string err_msg{"Could not find matching surface for periodic "
|
||||
"boundary condition on surface "};
|
||||
err_msg += std::to_string(surf->id);
|
||||
fatal_error(err_msg);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -1454,17 +1626,24 @@ surface_to_hdf5(int surf_ind, hid_t group)
|
|||
}
|
||||
|
||||
extern "C" bool
|
||||
surface_periodic(int surf_ind1, int surf_ind2, double xyz[3], double uvw[3])
|
||||
surface_periodic(int surf_ind1, double xyz[3], double uvw[3])
|
||||
{
|
||||
// Hopefully this function has only been called on a pair of surfaces that
|
||||
// Hopefully this function has only been called for a pair of surfaces that
|
||||
// support periodic BCs (checking should have been done when reading the
|
||||
// geometry XML). Downcast the surfaces to the PeriodicSurface type so we
|
||||
// can call the periodic_translate method.
|
||||
Surface *surf1_gen = surfaces_c[surf_ind1];
|
||||
PeriodicSurface *surf1 = dynamic_cast<PeriodicSurface *>(surf1_gen);
|
||||
Surface *surf2_gen = surfaces_c[surf_ind2];
|
||||
Surface *surf2_gen = surfaces_c[surf1->i_periodic];
|
||||
PeriodicSurface *surf2 = dynamic_cast<PeriodicSurface *>(surf2_gen);
|
||||
|
||||
// Call the type-bound methods.
|
||||
return surf2->periodic_translate(surf1, xyz, uvw);
|
||||
}
|
||||
|
||||
extern "C" int
|
||||
surface_i_periodic(int surf_ind)
|
||||
{
|
||||
PeriodicSurface *surf = dynamic_cast<PeriodicSurface *>(surfaces_c[surf_ind]);
|
||||
return surf->i_periodic;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -12,68 +12,17 @@ module surface_header
|
|||
! construct closed volumes (cells)
|
||||
!===============================================================================
|
||||
|
||||
type, abstract :: Surface
|
||||
type :: Surface
|
||||
integer :: id ! Unique ID
|
||||
integer, allocatable :: &
|
||||
neighbor_pos(:), & ! List of cells on positive side
|
||||
neighbor_neg(:) ! List of cells on negative side
|
||||
integer :: bc ! Boundary condition
|
||||
integer :: i_periodic = NONE ! Index of corresponding periodic surface
|
||||
character(len=104) :: name = "" ! User-defined name
|
||||
end type Surface
|
||||
|
||||
!===============================================================================
|
||||
! SURFACECONTAINER allows us to store an array of different types of surfaces
|
||||
!===============================================================================
|
||||
|
||||
type :: SurfaceContainer
|
||||
class(Surface), allocatable :: obj
|
||||
end type SurfaceContainer
|
||||
|
||||
!===============================================================================
|
||||
! All the derived types below are extensions of the abstract Surface type. They
|
||||
! inherent the reflect() type-bound procedure and must implement normal()
|
||||
!===============================================================================
|
||||
|
||||
type, extends(Surface) :: SurfaceXPlane
|
||||
end type SurfaceXPlane
|
||||
|
||||
type, extends(Surface) :: SurfaceYPlane
|
||||
end type SurfaceYPlane
|
||||
|
||||
type, extends(Surface) :: SurfaceZPlane
|
||||
end type SurfaceZPlane
|
||||
|
||||
type, extends(Surface) :: SurfacePlane
|
||||
end type SurfacePlane
|
||||
|
||||
type, extends(Surface) :: SurfaceXCylinder
|
||||
end type SurfaceXCylinder
|
||||
|
||||
type, extends(Surface) :: SurfaceYCylinder
|
||||
end type SurfaceYCylinder
|
||||
|
||||
type, extends(Surface) :: SurfaceZCylinder
|
||||
end type SurfaceZCylinder
|
||||
|
||||
type, extends(Surface) :: SurfaceSphere
|
||||
end type SurfaceSphere
|
||||
|
||||
type, extends(Surface) :: SurfaceXCone
|
||||
end type SurfaceXCone
|
||||
|
||||
type, extends(Surface) :: SurfaceYCone
|
||||
end type SurfaceYCone
|
||||
|
||||
type, extends(Surface) :: SurfaceZCone
|
||||
end type SurfaceZCone
|
||||
|
||||
type, extends(Surface) :: SurfaceQuadric
|
||||
end type SurfaceQuadric
|
||||
|
||||
integer(C_INT32_T), bind(C) :: n_surfaces ! # of surfaces
|
||||
|
||||
type(SurfaceContainer), allocatable, target :: surfaces(:)
|
||||
type(Surface), allocatable, target :: surfaces(:)
|
||||
|
||||
! Dictionary that maps user IDs to indices in 'surfaces'
|
||||
type(DictIntInt) :: surface_dict
|
||||
|
|
|
|||
|
|
@ -105,7 +105,7 @@ contains
|
|||
integer, intent(in) :: bin
|
||||
character(MAX_LINE_LEN) :: label
|
||||
|
||||
label = "Surface " // to_str(surfaces(this % surfaces(bin)) % obj % id)
|
||||
label = "Surface " // to_str(surfaces(this % surfaces(bin)) % id)
|
||||
end function text_label_surface
|
||||
|
||||
end module tally_filter_surface
|
||||
|
|
|
|||
|
|
@ -6,7 +6,7 @@ module tracking
|
|||
use geometry_header, only: cells
|
||||
use geometry, only: find_cell, distance_to_boundary, cross_lattice, &
|
||||
check_cell_overlap, surface_reflect_c, &
|
||||
surface_periodic_c
|
||||
surface_periodic_c, surface_i_periodic_c
|
||||
use message_passing
|
||||
use mgxs_header
|
||||
use nuclide_header
|
||||
|
|
@ -298,7 +298,7 @@ contains
|
|||
class(Surface), pointer :: surf
|
||||
|
||||
i_surface = abs(p % surface)
|
||||
surf => surfaces(i_surface)%obj
|
||||
surf => surfaces(i_surface)
|
||||
if (verbosity >= 10 .or. trace) then
|
||||
call write_message(" Crossing surface " // trim(to_str(surf % id)))
|
||||
end if
|
||||
|
|
@ -412,14 +412,14 @@ contains
|
|||
p % coord(1) % xyz = xyz
|
||||
end if
|
||||
|
||||
rotational = surface_periodic_c(i_surface-1, surf % i_periodic-1, &
|
||||
p % coord(1) % xyz, p % coord(1) % uvw)
|
||||
rotational = surface_periodic_c(i_surface-1, p % coord(1) % xyz, &
|
||||
p % coord(1) % uvw)
|
||||
|
||||
! Reassign particle's surface
|
||||
if (rotational) then
|
||||
p % surface = surf % i_periodic
|
||||
p % surface = surface_i_periodic_c(i_surface-1) + 1
|
||||
else
|
||||
p % surface = sign(surf % i_periodic, p % surface)
|
||||
p % surface = sign(surface_i_periodic_c(i_surface-1) + 1, p % surface)
|
||||
end if
|
||||
|
||||
! Figure out what cell particle is in now
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue