mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Use Position/Angle in lattice classes
This commit is contained in:
Paul Romano
parent
d9a7db1677
commit
f7649a61dc
4 changed files with 93 additions and 91 deletions
|
|
@ -12,7 +12,7 @@ struct Position {
|
|||
|
||||
Position() = default;
|
||||
Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { };
|
||||
Position(double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
|
||||
Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
|
||||
|
||||
Position& operator+=(Position);
|
||||
Position& operator+=(double);
|
||||
|
|
@ -59,7 +59,7 @@ struct Angle : Position {
|
|||
double& w() { return z; }
|
||||
Angle() = default;
|
||||
Angle(double u, double v, double w) : Position{u, v, w} { };
|
||||
Angle(double uvw[]) : Position{uvw} { };
|
||||
Angle(const double uvw[]) : Position{uvw} { };
|
||||
Angle(Position r) : Position{r} { };
|
||||
};
|
||||
|
||||
|
|
|
|||
|
|
@ -11,6 +11,8 @@
|
|||
#include <vector>
|
||||
#include <complex.h>
|
||||
|
||||
#include "geometry.h"
|
||||
|
||||
|
||||
namespace openmc {
|
||||
|
||||
|
|
@ -181,5 +183,12 @@ write_dataset(hid_t obj_id, const char* name, const std::array<T, N>& buffer)
|
|||
write_dataset(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data(), false);
|
||||
}
|
||||
|
||||
template<> inline void
|
||||
write_dataset<Position>(hid_t obj_id, const char* name, Position r)
|
||||
{
|
||||
std::array<double, 3> buffer {r.x, r.y, r.z};
|
||||
write_dataset(obj_id, name, buffer);
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
#endif //HDF5_INTERFACE_H
|
||||
|
|
|
|||
130
src/lattice.cpp
130
src/lattice.cpp
|
|
@ -223,15 +223,14 @@ RectLattice::are_valid_indices(const int i_xyz[3]) const
|
|||
//==============================================================================
|
||||
|
||||
std::pair<double, std::array<int, 3>>
|
||||
RectLattice::distance(const double xyz[3], const double uvw[3],
|
||||
const int i_xyz[3]) const
|
||||
RectLattice::distance(Position r, Angle a, const int i_xyz[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]};
|
||||
double x = r.x;
|
||||
double y = r.y;
|
||||
double z = r.z;
|
||||
double u = a.u();
|
||||
double v = a.v();
|
||||
|
||||
// Determine the oncoming edge.
|
||||
double x0 {copysign(0.5 * pitch[0], u)};
|
||||
|
|
@ -264,7 +263,7 @@ RectLattice::distance(const double xyz[3], const double uvw[3],
|
|||
|
||||
// Top and bottom sides
|
||||
if (is_3d) {
|
||||
double w {uvw[2]};
|
||||
double w {a.w()};
|
||||
double z0 {copysign(0.5 * pitch[2], w)};
|
||||
if ((std::abs(z - z0) > FP_PRECISION) && w != 0) {
|
||||
double this_d = (z0 - z) / w;
|
||||
|
|
@ -285,13 +284,13 @@ RectLattice::distance(const double xyz[3], const double uvw[3],
|
|||
//==============================================================================
|
||||
|
||||
std::array<int, 3>
|
||||
RectLattice::get_indices(const double xyz[3]) const
|
||||
RectLattice::get_indices(Position r) const
|
||||
{
|
||||
int ix {static_cast<int>(std::ceil((xyz[0] - lower_left[0]) / pitch[0]))-1};
|
||||
int iy {static_cast<int>(std::ceil((xyz[1] - lower_left[1]) / pitch[1]))-1};
|
||||
int ix {static_cast<int>(std::ceil((r.x - lower_left.x) / pitch.x))-1};
|
||||
int iy {static_cast<int>(std::ceil((r.y - lower_left.y) / pitch.y))-1};
|
||||
int iz;
|
||||
if (is_3d) {
|
||||
iz = static_cast<int>(std::ceil((xyz[2] - lower_left[2]) / pitch[2]))-1;
|
||||
iz = static_cast<int>(std::ceil((r.z - lower_left.z) / pitch.z))-1;
|
||||
} else {
|
||||
iz = 0;
|
||||
}
|
||||
|
|
@ -300,18 +299,15 @@ RectLattice::get_indices(const double xyz[3]) const
|
|||
|
||||
//==============================================================================
|
||||
|
||||
std::array<double, 3>
|
||||
RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
|
||||
Position
|
||||
RectLattice::get_local_position(Position r, const int i_xyz[3]) const
|
||||
{
|
||||
std::array<double, 3> local_xyz;
|
||||
local_xyz[0] = global_xyz[0] - (lower_left[0] + (i_xyz[0] + 0.5)*pitch[0]);
|
||||
local_xyz[1] = global_xyz[1] - (lower_left[1] + (i_xyz[1] + 0.5)*pitch[1]);
|
||||
r.x -= (lower_left.x + (i_xyz[0] + 0.5)*pitch.x);
|
||||
r.y -= (lower_left.y + (i_xyz[1] + 0.5)*pitch.y);
|
||||
if (is_3d) {
|
||||
local_xyz[2] = global_xyz[2] - (lower_left[2] + (i_xyz[2] + 0.5)*pitch[2]);
|
||||
} else {
|
||||
local_xyz[2] = global_xyz[2];
|
||||
r.z -= (lower_left.z + (i_xyz[2] + 0.5)*pitch.z);
|
||||
}
|
||||
return local_xyz;
|
||||
return r;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -583,12 +579,11 @@ HexLattice::are_valid_indices(const int i_xyz[3]) const
|
|||
//==============================================================================
|
||||
|
||||
std::pair<double, std::array<int, 3>>
|
||||
HexLattice::distance(const double xyz[3], const double uvw[3],
|
||||
const int i_xyz[3]) const
|
||||
HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const
|
||||
{
|
||||
// Compute the direction on the hexagonal basis.
|
||||
double beta_dir = uvw[0] * std::sqrt(3.0) / 2.0 + uvw[1] / 2.0;
|
||||
double gamma_dir = uvw[0] * std::sqrt(3.0) / 2.0 - uvw[1] / 2.0;
|
||||
double beta_dir = a.u() * std::sqrt(3.0) / 2.0 + a.v() / 2.0;
|
||||
double gamma_dir = a.u() * std::sqrt(3.0) / 2.0 - a.v() / 2.0;
|
||||
|
||||
// Note that hexagonal lattice distance calculations are performed
|
||||
// using the particle's coordinates relative to the neighbor lattice
|
||||
|
|
@ -600,15 +595,15 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
|
|||
double d {INFTY};
|
||||
std::array<int, 3> lattice_trans;
|
||||
double edge = -copysign(0.5*pitch[0], beta_dir); // Oncoming edge
|
||||
std::array<double, 3> xyz_t;
|
||||
Position r_t;
|
||||
if (beta_dir > 0) {
|
||||
const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1], i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
} else {
|
||||
const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1], i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
}
|
||||
double beta = xyz_t[0] * std::sqrt(3.0) / 2.0 + xyz_t[1] / 2.0;
|
||||
double beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0;
|
||||
if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) {
|
||||
d = (edge - beta) / beta_dir;
|
||||
if (beta_dir > 0) {
|
||||
|
|
@ -622,12 +617,12 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
|
|||
edge = -copysign(0.5*pitch[0], gamma_dir);
|
||||
if (gamma_dir > 0) {
|
||||
const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
} else {
|
||||
const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
}
|
||||
double gamma = xyz_t[0] * std::sqrt(3.0) / 2.0 - xyz_t[1] / 2.0;
|
||||
double gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0;
|
||||
if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) {
|
||||
double this_d = (edge - gamma) / gamma_dir;
|
||||
if (this_d < d) {
|
||||
|
|
@ -641,18 +636,18 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
|
|||
}
|
||||
|
||||
// Upper and lower sides.
|
||||
edge = -copysign(0.5*pitch[0], uvw[1]);
|
||||
if (uvw[1] > 0) {
|
||||
edge = -copysign(0.5*pitch[0], a.v());
|
||||
if (a.v() > 0) {
|
||||
const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
} else {
|
||||
const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]};
|
||||
xyz_t = get_local_xyz(xyz, i_xyz_t);
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
}
|
||||
if ((std::abs(xyz_t[1] - edge) > FP_PRECISION) && uvw[1] != 0) {
|
||||
double this_d = (edge - xyz_t[1]) / uvw[1];
|
||||
if ((std::abs(r_t.y - edge) > FP_PRECISION) && a.v() != 0) {
|
||||
double this_d = (edge - r_t.y) / a.v();
|
||||
if (this_d < d) {
|
||||
if (uvw[1] > 0) {
|
||||
if (a.v() > 0) {
|
||||
lattice_trans = {0, 1, 0};
|
||||
} else {
|
||||
lattice_trans = {0, -1, 0};
|
||||
|
|
@ -663,8 +658,8 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
|
|||
|
||||
// Top and bottom sides
|
||||
if (is_3d) {
|
||||
double z {xyz[2]};
|
||||
double w {uvw[2]};
|
||||
double z = r.z;
|
||||
double w = a.w();
|
||||
double z0 {copysign(0.5 * pitch[1], w)};
|
||||
if ((std::abs(z - z0) > FP_PRECISION) && w != 0) {
|
||||
double this_d = (z0 - z) / w;
|
||||
|
|
@ -686,24 +681,24 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
|
|||
//==============================================================================
|
||||
|
||||
std::array<int, 3>
|
||||
HexLattice::get_indices(const double xyz[3]) const
|
||||
HexLattice::get_indices(Position r) 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];}
|
||||
Position r_o {r.x - center.x, r.y - center.y, r.z};
|
||||
if (is_3d) {r_o.z -= center.z;}
|
||||
|
||||
// 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))-1;
|
||||
out[2] = static_cast<int>(std::ceil(r_o.z / pitch[1] + 0.5 * n_axial))-1;
|
||||
} else {
|
||||
out[2] = 0;
|
||||
}
|
||||
|
||||
// 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]
|
||||
double alpha = r_o.y - r_o.x / std::sqrt(3.0);
|
||||
out[0] = static_cast<int>(std::floor(r_o.x
|
||||
/ (0.5*std::sqrt(3.0) * pitch[0])));
|
||||
out[1] = static_cast<int>(std::floor(alpha / pitch[0]));
|
||||
|
||||
|
|
@ -725,8 +720,8 @@ HexLattice::get_indices(const double xyz[3]) const
|
|||
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];
|
||||
Position r_t = get_local_position(r, i_xyz);
|
||||
double d = r_t.x*r_t.x + r_t.y*r_t.y;
|
||||
if (d < d_min) {
|
||||
d_min = d;
|
||||
k_min = k;
|
||||
|
|
@ -751,26 +746,19 @@ HexLattice::get_indices(const double xyz[3]) const
|
|||
|
||||
//==============================================================================
|
||||
|
||||
std::array<double, 3>
|
||||
HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
|
||||
Position
|
||||
HexLattice::get_local_position(Position r, const int i_xyz[3]) const
|
||||
{
|
||||
std::array<double, 3> local_xyz;
|
||||
|
||||
// x_l = x_g - (center + pitch_x*cos(30)*index_x)
|
||||
local_xyz[0] = global_xyz[0] - (center[0]
|
||||
+ std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]);
|
||||
r.x -= (center.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]);
|
||||
// y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y)
|
||||
local_xyz[1] = global_xyz[1] - (center[1]
|
||||
+ (i_xyz[1] - n_rings + 1) * pitch[0]
|
||||
+ (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0);
|
||||
r.y -= (center.y + (i_xyz[1] - n_rings + 1) * pitch[0]
|
||||
+ (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0);
|
||||
if (is_3d) {
|
||||
local_xyz[2] = global_xyz[2] - center[2]
|
||||
+ (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1];
|
||||
} else {
|
||||
local_xyz[2] = global_xyz[2];
|
||||
r.z -= center.z - (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1];
|
||||
}
|
||||
|
||||
return local_xyz;
|
||||
return r;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -908,7 +896,9 @@ extern "C" {
|
|||
void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3],
|
||||
const int i_xyz[3], double *d, int lattice_trans[3])
|
||||
{
|
||||
std::pair<double, std::array<int, 3>> ld {lat->distance(xyz, uvw, i_xyz)};
|
||||
Position r {xyz};
|
||||
Angle a {uvw};
|
||||
std::pair<double, std::array<int, 3>> ld {lat->distance(r, a, i_xyz)};
|
||||
*d = ld.first;
|
||||
lattice_trans[0] = ld.second[0];
|
||||
lattice_trans[1] = ld.second[1];
|
||||
|
|
@ -917,7 +907,8 @@ extern "C" {
|
|||
|
||||
void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3])
|
||||
{
|
||||
std::array<int, 3> inds = lat->get_indices(xyz);
|
||||
Position r {xyz};
|
||||
std::array<int, 3> inds = lat->get_indices(r);
|
||||
i_xyz[0] = inds[0];
|
||||
i_xyz[1] = inds[1];
|
||||
i_xyz[2] = inds[2];
|
||||
|
|
@ -926,10 +917,11 @@ extern "C" {
|
|||
void lattice_get_local_xyz(Lattice *lat, const double global_xyz[3],
|
||||
const int i_xyz[3], double local_xyz[3])
|
||||
{
|
||||
std::array<double, 3> xyz = lat->get_local_xyz(global_xyz, i_xyz);
|
||||
local_xyz[0] = xyz[0];
|
||||
local_xyz[1] = xyz[1];
|
||||
local_xyz[2] = xyz[2];
|
||||
Position global {global_xyz};
|
||||
Position local = lat->get_local_position(global, i_xyz);
|
||||
local_xyz[0] = local.x;
|
||||
local_xyz[1] = local.y;
|
||||
local_xyz[2] = local.z;
|
||||
}
|
||||
|
||||
int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3])
|
||||
|
|
|
|||
|
|
@ -8,6 +8,7 @@
|
|||
#include <vector>
|
||||
|
||||
#include "constants.h"
|
||||
#include "geometry.h"
|
||||
#include "hdf5.h"
|
||||
#include "pugixml.hpp"
|
||||
|
||||
|
|
@ -75,26 +76,26 @@ public:
|
|||
virtual bool are_valid_indices(const int i_xyz[3]) const = 0;
|
||||
|
||||
//! \brief Find the next lattice surface crossing
|
||||
//! @param xyz[3] A 3D Cartesian coordinate.
|
||||
//! @param uvw[3] A 3D Cartesian direction.
|
||||
//! @param r A 3D Cartesian coordinate.
|
||||
//! @param a 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
|
||||
distance(Position r, Angle a, const int i_xyz[3]) const
|
||||
= 0;
|
||||
|
||||
//! \brief Find the lattice tile indices for a given point.
|
||||
//! @param xyz[3] A 3D Cartesian coordinate.
|
||||
//! @param r 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 std::array<int, 3> get_indices(Position r) const = 0;
|
||||
|
||||
//! \brief Get coordinates local to a lattice tile.
|
||||
//! @param global_xyz[3] A 3D Cartesian coordinate.
|
||||
//! @param r A 3D Cartesian coordinate.
|
||||
//! @param i_xyz[3] The indices for a lattice tile.
|
||||
//! @return Local 3D Cartesian coordinates.
|
||||
virtual std::array<double, 3>
|
||||
get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const = 0;
|
||||
virtual Position
|
||||
get_local_position(Position r, const int i_xyz[3]) const = 0;
|
||||
|
||||
//! \brief Check flattened lattice index.
|
||||
//! @param indx The index for a lattice tile.
|
||||
|
|
@ -193,12 +194,12 @@ public:
|
|||
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;
|
||||
distance(Position r, Angle a, const int i_xyz[3]) const;
|
||||
|
||||
std::array<int, 3> get_indices(const double xyz[3]) const;
|
||||
std::array<int, 3> get_indices(Position r) const;
|
||||
|
||||
std::array<double, 3>
|
||||
get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const;
|
||||
Position
|
||||
get_local_position(Position r, const int i_xyz[3]) const;
|
||||
|
||||
int32_t& offset(int map, const int i_xyz[3]);
|
||||
|
||||
|
|
@ -207,9 +208,9 @@ public:
|
|||
void to_hdf5_inner(hid_t group_id) const;
|
||||
|
||||
private:
|
||||
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
|
||||
std::array<int, 3> n_cells; //!< Number of cells along each axis
|
||||
Position lower_left; //!< Global lower-left corner of the lattice
|
||||
Position pitch; //!< Lattice tile width along each axis
|
||||
|
||||
// Convenience aliases
|
||||
int &nx {n_cells[0]};
|
||||
|
|
@ -233,12 +234,12 @@ public:
|
|||
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;
|
||||
distance(Position r, Angle a, const int i_xyz[3]) const;
|
||||
|
||||
std::array<int, 3> get_indices(const double xyz[3]) const;
|
||||
std::array<int, 3> get_indices(Position r) const;
|
||||
|
||||
std::array<double, 3>
|
||||
get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const;
|
||||
Position
|
||||
get_local_position(Position r, const int i_xyz[3]) const;
|
||||
|
||||
bool is_valid_index(int indx) const;
|
||||
|
||||
|
|
@ -251,7 +252,7 @@ public:
|
|||
private:
|
||||
int n_rings; //!< Number of radial tile positions
|
||||
int n_axial; //!< Number of axial tile positions
|
||||
std::array<double, 3> center; //!< Global center of lattice
|
||||
Position center; //!< Global center of lattice
|
||||
std::array<double, 2> pitch; //!< Lattice tile width and height
|
||||
};
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue