mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Add C++ RectLattice distance_to_boundary
This commit is contained in:
Sterling Harper
parent
98af5ed6be
commit
11336f875f
4 changed files with 127 additions and 73 deletions
|
|
@ -422,70 +422,8 @@ contains
|
|||
LAT_TYPE: select type(lat)
|
||||
|
||||
type is (RectLattice)
|
||||
! copy local coordinates
|
||||
x = p % coord(j) % xyz(1)
|
||||
y = p % coord(j) % xyz(2)
|
||||
z = p % coord(j) % xyz(3)
|
||||
|
||||
! determine oncoming edge
|
||||
x0 = sign(lat % pitch(1) * HALF, u)
|
||||
y0 = sign(lat % pitch(2) * HALF, v)
|
||||
|
||||
! left and right sides
|
||||
if (abs(x - x0) < FP_PRECISION) then
|
||||
d = INFINITY
|
||||
elseif (u == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
d = (x0 - x)/u
|
||||
end if
|
||||
|
||||
d_lat = d
|
||||
if (u > 0) then
|
||||
level_lat_trans(:) = [1, 0, 0]
|
||||
else
|
||||
level_lat_trans(:) = [-1, 0, 0]
|
||||
end if
|
||||
|
||||
! front and back sides
|
||||
if (abs(y - y0) < FP_PRECISION) then
|
||||
d = INFINITY
|
||||
elseif (v == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
d = (y0 - y)/v
|
||||
end if
|
||||
|
||||
if (d < d_lat) then
|
||||
d_lat = d
|
||||
if (v > 0) then
|
||||
level_lat_trans(:) = [0, 1, 0]
|
||||
else
|
||||
level_lat_trans(:) = [0, -1, 0]
|
||||
end if
|
||||
end if
|
||||
|
||||
if (lat % is_3d) then
|
||||
z0 = sign(lat % pitch(3) * HALF, w)
|
||||
|
||||
! top and bottom sides
|
||||
if (abs(z - z0) < FP_PRECISION) then
|
||||
d = INFINITY
|
||||
elseif (w == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
d = (z0 - z)/w
|
||||
end if
|
||||
|
||||
if (d < d_lat) then
|
||||
d_lat = d
|
||||
if (w > 0) then
|
||||
level_lat_trans(:) = [0, 0, 1]
|
||||
else
|
||||
level_lat_trans(:) = [0, 0, -1]
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
call lat % distance(p % coord(j) % xyz, p % coord(j) % uvw, &
|
||||
d_lat, level_lat_trans)
|
||||
|
||||
type is (HexLattice) LAT_TYPE
|
||||
! Copy local coordinates.
|
||||
|
|
|
|||
|
|
@ -95,6 +95,17 @@ module geometry_header
|
|||
integer(C_INT32_T) :: id
|
||||
end function lattice_id_c
|
||||
|
||||
subroutine lattice_distance_c(lat_ptr, xyz, uvw, d, lattice_trans) &
|
||||
bind(C, name='lattice_distance')
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
type(C_PTR), intent(in), value :: lat_ptr
|
||||
real(C_DOUBLE), intent(in) :: xyz(3)
|
||||
real(C_DOUBLE), intent(in) :: uvw(3)
|
||||
real(C_DOUBLE), intent(out) :: d
|
||||
integer(C_INT), intent(out) :: lattice_trans(3)
|
||||
end subroutine lattice_distance_c
|
||||
|
||||
subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5')
|
||||
use ISO_C_BINDING
|
||||
use hdf5
|
||||
|
|
@ -134,6 +145,7 @@ module geometry_header
|
|||
contains
|
||||
|
||||
procedure :: id => lattice_id
|
||||
procedure :: distance => lattice_distance
|
||||
procedure :: to_hdf5 => lattice_to_hdf5
|
||||
|
||||
procedure(lattice_are_valid_indices_), deferred :: are_valid_indices
|
||||
|
|
@ -281,6 +293,15 @@ contains
|
|||
id = lattice_id_c(this % ptr)
|
||||
end function lattice_id
|
||||
|
||||
subroutine lattice_distance(this, xyz, uvw, d, lattice_trans)
|
||||
class(Lattice), intent(in) :: this
|
||||
real(C_DOUBLE), intent(in) :: xyz(3)
|
||||
real(C_DOUBLE), intent(in) :: uvw(3)
|
||||
real(C_DOUBLE), intent(out) :: d
|
||||
integer(C_INT), intent(out) :: lattice_trans(3)
|
||||
call lattice_distance_c(this % ptr, xyz, uvw, d, lattice_trans)
|
||||
end subroutine lattice_distance
|
||||
|
||||
subroutine lattice_to_hdf5(this, group)
|
||||
class(Lattice), intent(in) :: this
|
||||
integer(HID_T), intent(in) :: group
|
||||
|
|
|
|||
|
|
@ -7,6 +7,9 @@
|
|||
#include "hdf5_interface.h"
|
||||
#include "xml_interface.h"
|
||||
|
||||
//TODO: this is only inlcuded for constants that should be moved elsewhere
|
||||
//#include "surface.h"
|
||||
|
||||
//TODO: remove this include
|
||||
#include <iostream>
|
||||
|
||||
|
|
@ -81,9 +84,9 @@ RectLattice::RectLattice(pugi::xml_node lat_node)
|
|||
fatal_error("Number of entries on <lower_left> must be the same as the "
|
||||
"number of entries on <dimension>.");
|
||||
}
|
||||
lower_left[0] = stoi(ll_words[0]);
|
||||
lower_left[1] = stoi(ll_words[1]);
|
||||
if (is_3d) {lower_left[2] = stoi(ll_words[2]);}
|
||||
lower_left[0] = stod(ll_words[0]);
|
||||
lower_left[1] = stod(ll_words[1]);
|
||||
if (is_3d) {lower_left[2] = stod(ll_words[2]);}
|
||||
|
||||
// Read the lattice pitches.
|
||||
std::string pitch_str{get_node_value(lat_node, "pitch")};
|
||||
|
|
@ -92,9 +95,9 @@ RectLattice::RectLattice(pugi::xml_node lat_node)
|
|||
fatal_error("Number of entries on <pitch> must be the same as the "
|
||||
"number of entries on <dimension>.");
|
||||
}
|
||||
pitch[0] = stoi(pitch_words[0]);
|
||||
pitch[1] = stoi(pitch_words[1]);
|
||||
if (is_3d) {pitch[2] = stoi(pitch_words[2]);}
|
||||
pitch[0] = stod(pitch_words[0]);
|
||||
pitch[1] = stod(pitch_words[1]);
|
||||
if (is_3d) {pitch[2] = stod(pitch_words[2]);}
|
||||
|
||||
// Read the universes and make sure the correct number was specified.
|
||||
int nx = n_cells[0];
|
||||
|
|
@ -129,6 +132,65 @@ RectLattice::RectLattice(pugi::xml_node lat_node)
|
|||
}
|
||||
}
|
||||
|
||||
std::pair<double, std::array<int, 3>>
|
||||
RectLattice::distance(const double xyz[3], const double uvw[3]) const
|
||||
{
|
||||
// Get short aliases to the coordinates.
|
||||
double x {xyz[0]};
|
||||
double y {xyz[1]};
|
||||
double z {xyz[2]};
|
||||
double u {uvw[0]};
|
||||
double v {uvw[1]};
|
||||
|
||||
// Determine the oncoming edge.
|
||||
double x0 {copysign(0.5 * pitch[0], u)};
|
||||
double y0 {copysign(0.5 * pitch[1], v)};
|
||||
|
||||
// Left and right sides
|
||||
double d {INFTY};
|
||||
std::array<int, 3> lattice_trans;
|
||||
if ((std::abs(x - x0) > FP_PRECISION) && u != 0) {
|
||||
d = (x0 - x) / u;
|
||||
if (u > 0) {
|
||||
lattice_trans = {1, 0, 0};
|
||||
} else {
|
||||
lattice_trans = {-1, 0, 0};
|
||||
}
|
||||
}
|
||||
|
||||
// Front and back sides
|
||||
if ((std::abs(y - y0) > FP_PRECISION) && v != 0) {
|
||||
double this_d = (y0 - y) / v;
|
||||
if (this_d < d) {
|
||||
d = this_d;
|
||||
if (v > 0) {
|
||||
lattice_trans = {0, 1, 0};
|
||||
} else {
|
||||
lattice_trans = {0, -1, 0};
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Top and bottom sides
|
||||
if (is_3d) {
|
||||
double w {uvw[2]};
|
||||
double z0 {copysign(0.5 * pitch[2], w)};
|
||||
if ((std::abs(z - z0) > FP_PRECISION) && w != 0) {
|
||||
double this_d = (z0 - z) / w;
|
||||
if (this_d < d) {
|
||||
d = this_d;
|
||||
if (w > 0) {
|
||||
lattice_trans = {0, 0, 1};
|
||||
} else {
|
||||
lattice_trans = {0, 0, -1};
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return {d, lattice_trans};
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// HexLattice implementation
|
||||
//==============================================================================
|
||||
|
|
@ -138,6 +200,11 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
|
|||
{
|
||||
}
|
||||
|
||||
std::pair<double, std::array<int, 3>>
|
||||
HexLattice::distance(const double xyz[3], const double uvw[3]) const
|
||||
{
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
||||
extern "C" void
|
||||
|
|
@ -173,6 +240,16 @@ extern "C" {
|
|||
|
||||
int32_t lattice_id(Lattice *lat) {return lat->id;}
|
||||
|
||||
void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3],
|
||||
double *d, int lattice_trans[3])
|
||||
{
|
||||
std::pair<double, std::array<int, 3>> ld {lat->distance(xyz, uvw)};
|
||||
*d = ld.first;
|
||||
lattice_trans[0] = ld.second[0];
|
||||
lattice_trans[1] = ld.second[1];
|
||||
lattice_trans[2] = ld.second[2];
|
||||
}
|
||||
|
||||
void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -3,6 +3,7 @@
|
|||
|
||||
#include <array>
|
||||
#include <cstdint>
|
||||
#include <limits> // For numeric_limits
|
||||
#include <map>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
|
@ -13,6 +14,14 @@
|
|||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// Constants that should eventually be moved out of this file
|
||||
//==============================================================================
|
||||
|
||||
extern "C" double FP_PRECISION;
|
||||
constexpr double INFTY{std::numeric_limits<double>::max()};
|
||||
constexpr int C_NONE {-1};
|
||||
|
||||
//==============================================================================
|
||||
// Constants
|
||||
//==============================================================================
|
||||
|
|
@ -51,6 +60,9 @@ public:
|
|||
|
||||
//virtual bool are_valid_indices(const int i_xyz[3]) const = 0;
|
||||
|
||||
virtual std::pair<double, std::array<int, 3>>
|
||||
distance(const double xyz[3], const double uvw[3]) const = 0;
|
||||
|
||||
//virtual void get_indices(const double global_xyz[3], int i_xyz[3]) const = 0;
|
||||
|
||||
//virtual void get_local_xyz(const double global_xyz[3], const int i_xyz[3],
|
||||
|
|
@ -74,10 +86,13 @@ public:
|
|||
|
||||
virtual ~RectLattice() {}
|
||||
|
||||
std::pair<double, std::array<int, 3>>
|
||||
distance(const double xyz[3], const double uvw[3]) const;
|
||||
|
||||
protected:
|
||||
std::array<int, 3> n_cells; //! Number of cells along each axis
|
||||
std::array<int, 3> lower_left; //! Global lower-left corner of the lattice
|
||||
std::array<int, 3> pitch; //! Lattice tile width along each axis
|
||||
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
|
||||
std::array<double, 3> pitch; //! Lattice tile width along each axis
|
||||
};
|
||||
|
||||
class HexLattice : public Lattice
|
||||
|
|
@ -86,6 +101,9 @@ public:
|
|||
explicit HexLattice(pugi::xml_node lat_node);
|
||||
|
||||
virtual ~HexLattice() {}
|
||||
|
||||
std::pair<double, std::array<int, 3>>
|
||||
distance(const double xyz[3], const double uvw[3]) const;
|
||||
};
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue