mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
lattice translation with S.Harper basis
This commit is contained in:
dr.yuri
parent
7059b56e79
commit
4ef3e30755
1 changed files with 127 additions and 124 deletions
251
src/lattice.cpp
251
src/lattice.cpp
|
|
@ -694,134 +694,137 @@ std::pair<double, std::array<int, 3>>
|
|||
HexLattice::distance(Position r, Direction u, const std::array<int, 3>& i_xyz)
|
||||
const
|
||||
{
|
||||
double cosa,sina,u_xy;
|
||||
// Index + 1 of translation direction
|
||||
int i_trans {1};
|
||||
xt :: xarray<int> translation_matrix;
|
||||
if (hextype==0){
|
||||
cosa = std::sqrt(3.0) / 2.0;
|
||||
sina = 0.5;
|
||||
u_xy = u.y;
|
||||
//Translation matrix 4x3 {{beta},{gamma},{oy},{oz}} for OY case
|
||||
translation_matrix = {{1, 0, 0},
|
||||
{1, -1, 0},
|
||||
{0, 1, 0},
|
||||
{0, 0, 1}};
|
||||
} else {
|
||||
cosa = 0.5;
|
||||
sina = std::sqrt(3.0) / 2.0;
|
||||
u_xy = u.x;
|
||||
//Translation matrix 4x3 {{beta},{gamma},{ox},{oz}} for OX case
|
||||
translation_matrix = {{0, 1, 0},
|
||||
{1, -1, 0},
|
||||
{1, 0, 0},
|
||||
{0, 0, 1}};
|
||||
}
|
||||
double beta_dir = u.x * cosa + u.y * sina;
|
||||
double gamma_dir = u.x * cosa - u.y * sina;
|
||||
//Short description
|
||||
//OY - orientation:
|
||||
// basis0 = (1, 0)
|
||||
// basis1 = (-1/sqrt(3), 1) = +120 degrees from basis0
|
||||
// beta = (sqrt(3)/2, 1/2) = +30 degrees from basis0
|
||||
// gamma = (sqrt(3)/2, -1/2) = -60 degrees from beta
|
||||
// delta = (0, 1) = +60 degrees from beta
|
||||
//OX - orientation:
|
||||
// basis0 = (1/sqrt(3), -1)
|
||||
// basis1 = (0, 1) = +120 degrees from basis0
|
||||
// beta = (1, 0) = +30 degrees from basis0
|
||||
// gamma = (1/2, -sqrt(3)/2) = -60 degrees from beta
|
||||
// delta = (1/2, sqrt(3)/2) = +60 degrees from beta
|
||||
//OZ be considered separetly
|
||||
double beta_dir;
|
||||
double gamma_dir;
|
||||
double delta_dir;
|
||||
if (hextype==0){
|
||||
beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0;
|
||||
gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0;
|
||||
delta_dir = u.y;
|
||||
}
|
||||
else{
|
||||
beta_dir = u.x;
|
||||
gamma_dir = u.x / 2.0 - u.y * std::sqrt(3.0) / 2.0;
|
||||
delta_dir = u.x / 2.0 + u.y * std::sqrt(3.0) / 2.0;
|
||||
}
|
||||
|
||||
// Note that hexagonal lattice distance calculations are performed
|
||||
// using the particle's coordinates relative to the neighbor lattice
|
||||
// cells, not relative to the particle's current cell. This is done
|
||||
// because there is significant disagreement between neighboring cells
|
||||
// on where the lattice boundary is due to finite precision issues.
|
||||
// Note that hexagonal lattice distance calculations are performed
|
||||
// using the particle's coordinates relative to the neighbor lattice
|
||||
// cells, not relative to the particle's current cell. This is done
|
||||
// because there is significant disagreement between neighboring cells
|
||||
// on where the lattice boundary is due to finite precision issues.
|
||||
|
||||
//beta direction
|
||||
double d {INFTY};
|
||||
std::array<int, 3> lattice_trans;
|
||||
double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge
|
||||
Position r_t;
|
||||
if (beta_dir > 0) {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0]+translation_matrix(0,0), i_xyz[1]+translation_matrix(0,1), i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
} else {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0]-translation_matrix(0,0), i_xyz[1]-translation_matrix(0,1), i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
}
|
||||
double beta = r_t.x * cosa + r_t.y * sina;
|
||||
if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) {
|
||||
d = (edge - beta) / beta_dir;
|
||||
if (beta_dir > 0) {
|
||||
i_trans = 1;
|
||||
} else {
|
||||
i_trans = -1;
|
||||
}
|
||||
}
|
||||
//gamma direction
|
||||
edge = -copysign(0.5*pitch_[0], gamma_dir);
|
||||
if (gamma_dir > 0) {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0]+translation_matrix(1,0), i_xyz[1]+translation_matrix(1,1), i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
} else {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0]-translation_matrix(1,0), i_xyz[1]-translation_matrix(1,1), i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
}
|
||||
double gamma = r_t.x * cosa - r_t.y * sina;
|
||||
if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) {
|
||||
double this_d = (edge - gamma) / gamma_dir;
|
||||
if (this_d < d) {
|
||||
if (gamma_dir > 0) {
|
||||
i_trans = 2;
|
||||
} else {
|
||||
i_trans = -2;
|
||||
}
|
||||
d = this_d;
|
||||
}
|
||||
}
|
||||
//y or x direction
|
||||
edge = -copysign(0.5*pitch_[0], u_xy);
|
||||
if (u_xy > 0) {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0]+translation_matrix(2,0), i_xyz[1]+translation_matrix(2,1), i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
} else {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0]-translation_matrix(2,0), i_xyz[1]-translation_matrix(2,1), i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
}
|
||||
double r_txy;
|
||||
if (hextype==0){
|
||||
r_txy = r_t.y;
|
||||
}
|
||||
else {
|
||||
r_txy = r_t.x;
|
||||
}
|
||||
if ((std::abs(r_txy - edge) > FP_PRECISION) && u_xy != 0) {
|
||||
double this_d = (edge - r_txy) / u_xy;
|
||||
if (this_d < d) {
|
||||
if (u_xy > 0) {
|
||||
i_trans = 3;
|
||||
} else {
|
||||
i_trans = -3;
|
||||
}
|
||||
d = this_d;
|
||||
}
|
||||
}
|
||||
// Top and bottom sides
|
||||
if (is_3d_) {
|
||||
double z = r.z;
|
||||
double z0 {copysign(0.5 * pitch_[1], u.z)};
|
||||
if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) {
|
||||
double this_d = (z0 - z) / u.z;
|
||||
if (this_d < d) {
|
||||
d = this_d;
|
||||
if (u.z > 0) {
|
||||
i_trans = 4;
|
||||
} else {
|
||||
i_trans = -4;
|
||||
}
|
||||
d = this_d;
|
||||
}
|
||||
}
|
||||
}
|
||||
//beta direction
|
||||
double d {INFTY};
|
||||
std::array<int, 3> lattice_trans;
|
||||
double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge
|
||||
double beta;
|
||||
Position r_t;
|
||||
if (beta_dir > 0) {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0]+1, i_xyz[1], i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
} else {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0]-1, i_xyz[1], i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
}
|
||||
if (hextype==0){
|
||||
beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0;}
|
||||
else{
|
||||
beta = r_t.x ;}
|
||||
if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) {
|
||||
d = (edge - beta) / beta_dir;
|
||||
if (beta_dir > 0) {
|
||||
lattice_trans = {1, 0, 0};
|
||||
} else {
|
||||
lattice_trans = {-1, 0, 0};
|
||||
}
|
||||
}
|
||||
|
||||
for (int i=0;i<3;i++){
|
||||
if (i_trans < 0){
|
||||
lattice_trans[i] = -translation_matrix(-(i_trans+1),i);
|
||||
}else{
|
||||
lattice_trans[i] = translation_matrix((i_trans-1),i);
|
||||
}
|
||||
}
|
||||
// gamma direction.
|
||||
edge = -copysign(0.5*pitch_[0], gamma_dir);
|
||||
if (gamma_dir > 0) {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
} else {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
}
|
||||
double gamma;
|
||||
if (hextype==0){
|
||||
gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0;}
|
||||
else{
|
||||
gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0;}
|
||||
if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) {
|
||||
double this_d = (edge - gamma) / gamma_dir;
|
||||
if (this_d < d) {
|
||||
if (gamma_dir > 0) {
|
||||
lattice_trans = {1, -1, 0};
|
||||
} else {
|
||||
lattice_trans = {-1, 1, 0};
|
||||
}
|
||||
d = this_d;
|
||||
}
|
||||
}
|
||||
|
||||
return {d, lattice_trans};
|
||||
// delta directions.
|
||||
edge = -copysign(0.5*pitch_[0], delta_dir);
|
||||
if (delta_dir > 0) {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0], i_xyz[1]+1, i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
} else {
|
||||
const std::array<int, 3> i_xyz_t {i_xyz[0], i_xyz[1]-1, i_xyz[2]};
|
||||
r_t = get_local_position(r, i_xyz_t);
|
||||
}
|
||||
double delta;
|
||||
if (hextype==0){
|
||||
delta = r_t.y ;}
|
||||
else{
|
||||
delta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0;}
|
||||
if ((std::abs(delta - edge) > FP_PRECISION) && delta_dir != 0) {
|
||||
double this_d = (edge - delta) / delta_dir;
|
||||
if (this_d < d) {
|
||||
if (delta_dir > 0) {
|
||||
lattice_trans = {0, 1, 0};
|
||||
} else {
|
||||
lattice_trans = {0, -1, 0};
|
||||
}
|
||||
d = this_d;
|
||||
}
|
||||
}
|
||||
|
||||
// Top and bottom sides
|
||||
if (is_3d_) {
|
||||
double z = r.z;
|
||||
double z0 {copysign(0.5 * pitch_[1], u.z)};
|
||||
if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) {
|
||||
double this_d = (z0 - z) / u.z;
|
||||
if (this_d < d) {
|
||||
d = this_d;
|
||||
if (u.z > 0) {
|
||||
lattice_trans = {0, 0, 1};
|
||||
} else {
|
||||
lattice_trans = {0, 0, -1};
|
||||
}
|
||||
d = this_d;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return {d, lattice_trans};
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue