crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Move lattice get_indices and are_valid to C++

Browse source
This commit is contained in:
Sterling Harper 2018-02-16 17:49:42 -05:00
parent ec14970caf
commit efc77b74fe
3 changed files with 178 additions and 147 deletions
Download patch file Download diff file Expand all files Collapse all files

177
src/geometry_header.F90
View file

@ -95,6 +95,15 @@ module geometry_header
integer(C_INT32_T) :: id
end function lattice_id_c
function lattice_are_valid_indices_c(lat_ptr, i_xyz) &
bind(C, name='lattice_are_valid_indices') result (is_valid)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: lat_ptr
integer(C_INT), intent(in) :: i_xyz(3)
logical(C_BOOL) :: is_valid
end function lattice_are_valid_indices_c
subroutine lattice_distance_c(lat_ptr, xyz, uvw, i_xyz, d, lattice_trans) &
bind(C, name='lattice_distance')
use ISO_C_BINDING
@ -107,6 +116,15 @@ module geometry_header
integer(C_INT), intent(out) :: lattice_trans(3)
end subroutine lattice_distance_c
subroutine lattice_get_indices_c(lat_ptr, xyz, i_xyz) &
bind(C, name='lattice_get_indices')
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: lat_ptr
real(C_DOUBLE), intent(in) :: xyz(3)
integer(C_INT), intent(out) :: i_xyz(3)
end subroutine lattice_get_indices_c
subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5')
use ISO_C_BINDING
use hdf5
@ -146,37 +164,16 @@ module geometry_header
contains
procedure :: id => lattice_id
procedure :: are_valid_indices => lattice_are_valid_indices
procedure :: distance => lattice_distance
procedure :: get_indices => lattice_get_indices
procedure :: to_hdf5 => lattice_to_hdf5
procedure(lattice_are_valid_indices_), deferred :: are_valid_indices
procedure(lattice_get_indices_), deferred :: get_indices
procedure(lattice_get_local_xyz_), deferred :: get_local_xyz
end type Lattice
abstract interface
!===============================================================================
! ARE_VALID_INDICES returns .true. if the given lattice indices fit within the
! bounds of the lattice. Returns false otherwise.
function lattice_are_valid_indices_(this, i_xyz) result(is_valid)
import Lattice
class(Lattice), intent(in) :: this
integer, intent(in) :: i_xyz(3)
logical :: is_valid
end function lattice_are_valid_indices_
!===============================================================================
! GET_INDICES returns the indices in a lattice for the given global xyz.
function lattice_get_indices_(this, global_xyz) result(i_xyz)
import Lattice
class(Lattice), intent(in) :: this
real(8), intent(in) :: global_xyz(3)
integer :: i_xyz(3)
end function lattice_get_indices_
!===============================================================================
! GET_LOCAL_XYZ returns the translated local version of the given global xyz.
@ -199,8 +196,6 @@ module geometry_header
contains
procedure :: are_valid_indices => valid_inds_rect
procedure :: get_indices => get_inds_rect
procedure :: get_local_xyz => get_local_rect
end type RectLattice
@ -215,8 +210,6 @@ module geometry_header
contains
procedure :: are_valid_indices => valid_inds_hex
procedure :: get_indices => get_inds_hex
procedure :: get_local_xyz => get_local_hex
end type HexLattice
@ -294,6 +287,13 @@ contains
id = lattice_id_c(this % ptr)
end function lattice_id
function lattice_are_valid_indices(this, i_xyz) result (is_valid)
class(Lattice), intent(in) :: this
integer(C_INT), intent(in) :: i_xyz(3)
logical(C_BOOL) :: is_valid
is_valid = lattice_are_valid_indices_c(this % ptr, i_xyz)
end function lattice_are_valid_indices
subroutine lattice_distance(this, xyz, uvw, i_xyz, d, lattice_trans)
class(Lattice), intent(in) :: this
real(C_DOUBLE), intent(in) :: xyz(3)
@ -304,128 +304,19 @@ contains
call lattice_distance_c(this % ptr, xyz, uvw, i_xyz, d, lattice_trans)
end subroutine lattice_distance
function lattice_get_indices(this, xyz) result(i_xyz)
class(Lattice), intent(in) :: this
real(C_DOUBLE), intent(in) :: xyz(3)
integer(C_INT) :: i_xyz(3)
call lattice_get_indices_c(this % ptr, xyz, i_xyz)
end function lattice_get_indices
subroutine lattice_to_hdf5(this, group)
class(Lattice), intent(in) :: this
integer(HID_T), intent(in) :: group
call lattice_to_hdf5_c(this % ptr, group)
end subroutine lattice_to_hdf5
!===============================================================================
function valid_inds_rect(this, i_xyz) result(is_valid)
class(RectLattice), intent(in) :: this
integer, intent(in) :: i_xyz(3)
logical :: is_valid
is_valid = all(i_xyz > 0 .and. i_xyz <= this % n_cells)
end function valid_inds_rect
!===============================================================================
function valid_inds_hex(this, i_xyz) result(is_valid)
class(HexLattice), intent(in) :: this
integer, intent(in) :: i_xyz(3)
logical :: is_valid
is_valid = (all(i_xyz > 0) .and. &
&i_xyz(1) < 2*this % n_rings .and. &
&i_xyz(2) < 2*this % n_rings .and. &
&i_xyz(1) + i_xyz(2) > this % n_rings .and. &
&i_xyz(1) + i_xyz(2) < 3*this % n_rings .and. &
&i_xyz(3) <= this % n_axial)
end function valid_inds_hex
!===============================================================================
function get_inds_rect(this, global_xyz) result(i_xyz)
class(RectLattice), intent(in) :: this
real(8), intent(in) :: global_xyz(3)
integer :: i_xyz(3)
real(8) :: xyz(3) ! global_xyz alias
xyz = global_xyz
i_xyz(1) = ceiling((xyz(1) - this % lower_left(1))/this % pitch(1))
i_xyz(2) = ceiling((xyz(2) - this % lower_left(2))/this % pitch(2))
if (this % is_3d) then
i_xyz(3) = ceiling((xyz(3) - this % lower_left(3))/this % pitch(3))
else
i_xyz(3) = 1
end if
end function get_inds_rect
!===============================================================================
function get_inds_hex(this, global_xyz) result(i_xyz)
class(HexLattice), intent(in) :: this
real(8), intent(in) :: global_xyz(3)
integer :: i_xyz(3)
real(8) :: xyz(3) ! global xyz relative to the center
real(8) :: alpha ! Skewed coord axis
real(8) :: xyz_t(3) ! Local xyz
real(8) :: d, d_min ! Squared distance from cell centers
integer :: i, j, k ! Iterators
integer :: k_min ! Minimum distance index
xyz(1) = global_xyz(1) - this % center(1)
xyz(2) = global_xyz(2) - this % center(2)
! Index z direction.
if (this % is_3d) then
xyz(3) = global_xyz(3) - this % center(3)
i_xyz(3) = ceiling(xyz(3)/this % pitch(2) + HALF*this % n_axial)
else
xyz(3) = global_xyz(3)
i_xyz(3) = 1
end if
! Convert coordinates into skewed bases. The (x, alpha) basis is used to
! find the index of the global coordinates to within 4 cells.
alpha = xyz(2) - xyz(1) / sqrt(THREE)
i_xyz(1) = floor(xyz(1) / (sqrt(THREE) / TWO * this % pitch(1)))
i_xyz(2) = floor(alpha / this % pitch(1))
! Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
! the array is offset so that the indices never go below 1).
i_xyz(1) = i_xyz(1) + this % n_rings
i_xyz(2) = i_xyz(2) + this % n_rings
! Calculate the (squared) distance between the particle and the centers of
! the four possible cells. Regular hexagonal tiles form a centroidal
! Voronoi tessellation so the global xyz should be in the hexagonal cell
! that it is closest to the center of. This method is used over a
! method that uses the remainders of the floor divisions above because it
! provides better finite precision performance. Squared distances are
! used becasue they are more computationally efficient than normal
! distances.
k = 1
d_min = INFINITY
do i = 0, 1
do j = 0, 1
xyz_t = this % get_local_xyz(global_xyz, i_xyz + [j, i, 0])
d = xyz_t(1)**2 + xyz_t(2)**2
if (d < d_min) then
d_min = d
k_min = k
end if
k = k + 1
end do
end do
! Select the minimum squared distance which corresponds to the cell the
! coordinates are in.
if (k_min == 2) then
i_xyz(1) = i_xyz(1) + 1
else if (k_min == 3) then
i_xyz(2) = i_xyz(2) + 1
else if (k_min == 4) then
i_xyz(1) = i_xyz(1) + 1
i_xyz(2) = i_xyz(2) + 1
end if
end function get_inds_hex
!===============================================================================
function get_local_rect(this, global_xyz, i_xyz) result(local_xyz)

121
src/lattice.cpp
View file

@ -131,6 +131,18 @@ RectLattice::RectLattice(pugi::xml_node lat_node)
}
}
//==============================================================================
bool
RectLattice::are_valid_indices(const int i_xyz[3]) const
{
return ( (i_xyz[0] > 0) && (i_xyz[0] <= n_cells[0])
&& (i_xyz[1] > 0) && (i_xyz[1] <= n_cells[1])
&& (i_xyz[2] > 0) && (i_xyz[2] <= n_cells[2]));
}
//==============================================================================
std::pair<double, std::array<int, 3>>
RectLattice::distance(const double xyz[3], const double uvw[3],
const int i_xyz[3]) const
@ -191,6 +203,22 @@ RectLattice::distance(const double xyz[3], const double uvw[3],
return {d, lattice_trans};
}
//==============================================================================
std::array<int, 3>
RectLattice::get_indices(const double xyz[3]) const
{
int ix {static_cast<int>(std::ceil((xyz[0] - lower_left[0]) / pitch[0]))};
int iy {static_cast<int>(std::ceil((xyz[1] - lower_left[1]) / pitch[1]))};
int iz;
if (is_3d) {
iz = static_cast<int>(std::ceil((xyz[2] - lower_left[2]) / pitch[2]));
} else {
iz = 1;
}
return {ix, iy, iz};
}
//==============================================================================
// HexLattice implementation
//==============================================================================
@ -229,7 +257,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
fatal_error("A hexagonal lattice with <n_axial> must have <pitch> "
"specified by 2 numbers.");
} else if (!is_3d && (pitch_words.size() != 1)) {
fatal_error("A hexagonal lattice without <n_axial> must have <center> "
fatal_error("A hexagonal lattice without <n_axial> must have <pitch> "
"specified by 1 number.");
}
pitch[0] = stod(pitch_words[0]);
@ -342,6 +370,18 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
}
}
//==============================================================================
bool
HexLattice::are_valid_indices(const int i_xyz[3]) const
{
return ((i_xyz[0] > 0) && (i_xyz[1] > 0) && (i_xyz[2] > 0)
&& (i_xyz[0] < 2*n_rings) && (i_xyz[1] < 2*n_rings)
&& (i_xyz[0] + i_xyz[1] > n_rings)
&& (i_xyz[0] + i_xyz[1] < 3*n_rings)
&& (i_xyz[2] <= n_axial));
}
std::pair<double, std::array<int, 3>>
HexLattice::distance(const double xyz[3], const double uvw[3],
const int i_xyz[3]) const
@ -443,6 +483,74 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
return {d, lattice_trans};
}
//==============================================================================
std::array<int, 3>
HexLattice::get_indices(const double xyz[3]) const
{
// Offset the xyz by the lattice center.
double xyz_o[3] {xyz[0] - center[0], xyz[1] - center[1], xyz[2]};
if (is_3d) {xyz_o[2] -= center[2];}
// Index the z direction.
std::array<int, 3> out;
if (is_3d) {
out[2] = static_cast<int>(std::ceil(xyz_o[2] / pitch[1] + 0.5 * n_axial));
} else {
out[2] = 1;
}
// Convert coordinates into skewed bases. The (x, alpha) basis is used to
// find the index of the global coordinates to within 4 cells.
double alpha = xyz_o[1] - xyz_o[0] / std::sqrt(3.0);
out[0] = static_cast<int>(std::floor(xyz_o[0]
/ (0.5*std::sqrt(3.0) * pitch[0])));
out[1] = static_cast<int>(std::floor(alpha / pitch[0]));
// Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
// the array is offset so that the indices never go below 1).
out[0] += n_rings;
out[1] += n_rings;
// Calculate the (squared) distance between the particle and the centers of
// the four possible cells. Regular hexagonal tiles form a Voronoi
// tessellation so the xyz should be in the hexagonal cell that it is closest
// to the center of. This method is used over a method that uses the
// remainders of the floor divisions above because it provides better finite
// precision performance. Squared distances are used becasue they are more
// computationally efficient than normal distances.
int k {1};
int k_min {1};
double d_min {INFTY};
for (int i = 0; i < 2; i++) {
for (int j = 0; j < 2; j++) {
int i_xyz[3] {out[0] + j, out[1] + i, 0};
std::array<double, 3> xyz_t = get_local_xyz(xyz, i_xyz);
double d = xyz_t[0]*xyz_t[0] + xyz_t[1]*xyz_t[1];
if (d < d_min) {
d_min = d;
k_min = k;
}
k++;
}
}
// Select the minimum squared distance which corresponds to the cell the
// coordinates are in.
if (k_min == 2) {
out[0] += 1;
} else if (k_min == 3) {
out[1] += 1;
} else if (k_min == 4) {
out[0] += 1;
out[1] += 1;
}
return out;
}
//==============================================================================
std::array<double, 3>
HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
{
@ -500,6 +608,9 @@ extern "C" {
int32_t lattice_id(Lattice *lat) {return lat->id;}
bool lattice_are_valid_indices(Lattice *lat, const int i_xyz[3])
{return lat->are_valid_indices(i_xyz);}
void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3],
const int i_xyz[3], double *d, int lattice_trans[3])
{
@ -510,6 +621,14 @@ extern "C" {
lattice_trans[2] = ld.second[2];
}
void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3])
{
std::array<int, 3> inds {lat->get_indices(xyz)};
i_xyz[0] = inds[0];
i_xyz[1] = inds[1];
i_xyz[2] = inds[2];
}
void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);}
}

27
src/lattice.h
View file

@ -58,13 +58,26 @@ public:
virtual ~Lattice() {}
//virtual bool are_valid_indices(const int i_xyz[3]) const = 0;
//! Check lattice indices.
//! @param i_xyz[3] The indices for a lattice tile.
//! @return true if the given indices fit within the lattice bounds. False
//! otherwise.
virtual bool are_valid_indices(const int i_xyz[3]) const = 0;
//! Find the next lattice surface crossing
//! @param xyz[3] A 3D Cartesian coordinate.
//! @param uvw[3] A 3D Cartesian direction.
//! @param i_xyz[3] The indices for a lattice tile.
//! @return The distance to the next crossing and an array indicating how the
//! lattice indices would change after crossing that boundary.
virtual std::pair<double, std::array<int, 3>>
distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const
= 0;
= 0;
//virtual void get_indices(const double global_xyz[3], int i_xyz[3]) const = 0;
//! Find the lattice tile indices for a given point.
//! @param xyz[3] A 3D Cartesian coordinate.
//! @return An array containing the indices of a lattice tile.
virtual std::array<int, 3> get_indices(const double xyz[3]) const = 0;
//virtual void get_local_xyz(const double global_xyz[3], const int i_xyz[3],
// double local_xyz[3]) const = 0;
@ -87,9 +100,13 @@ public:
virtual ~RectLattice() {}
bool are_valid_indices(const int i_xyz[3]) const;
std::pair<double, std::array<int, 3>>
distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const;
std::array<int, 3> get_indices(const double xyz[3]) const;
protected:
std::array<int, 3> n_cells; //! Number of cells along each axis
std::array<double, 3> lower_left; //! Global lower-left corner of the lattice
@ -103,9 +120,13 @@ public:
virtual ~HexLattice() {}
bool are_valid_indices(const int i_xyz[3]) const;
std::pair<double, std::array<int, 3>>
distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const;
std::array<int, 3> get_indices(const double xyz[3]) const;
std::array<double, 3>
get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const;