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

Settle on Position/Direction, declared in position.h

Browse source
This commit is contained in:
Paul Romano 2018-08-08 13:36:42 -05:00
parent f7649a61dc
commit 0fba0e287c
11 changed files with 275 additions and 257 deletions
Download patch file Download diff file Expand all files Collapse all files

2
CMakeLists.txt
View file

@ -386,7 +386,6 @@ add_library(libopenmc SHARED
src/cell.cpp
src/initialize.cpp
src/finalize.cpp
src/geometry.cpp
src/geometry_aux.cpp
src/hdf5_interface.cpp
src/lattice.cpp
@ -396,6 +395,7 @@ add_library(libopenmc SHARED
src/mgxs_interface.cpp
src/particle.cpp
src/plot.cpp
src/position.cpp
src/pugixml/pugixml_c.cpp
src/random_lcg.cpp
src/scattdata.cpp

26
src/cell.cpp
View file

@ -279,19 +279,19 @@ Cell::Cell(pugi::xml_node cell_node)
//==============================================================================
bool
Cell::contains(Position r, Angle a, int32_t on_surface) const
Cell::contains(Position r, Direction u, int32_t on_surface) const
{
if (simple) {
return contains_simple(r, a, on_surface);
return contains_simple(r, u, on_surface);
} else {
return contains_complex(r, a, on_surface);
return contains_complex(r, u, on_surface);
}
}
//==============================================================================
std::pair<double, int32_t>
Cell::distance(Position r, Angle a, int32_t on_surface) const
Cell::distance(Position r, Direction u, int32_t on_surface) const
{
double min_dist {INFTY};
int32_t i_surf {std::numeric_limits<int32_t>::max()};
@ -303,7 +303,7 @@ Cell::distance(Position r, Angle a, int32_t on_surface) const
// Calculate the distance to this surface.
// Note the off-by-one indexing
bool coincident {token == on_surface};
double d {surfaces_c[abs(token)-1]->distance(r, a, coincident)};
double d {surfaces_c[abs(token)-1]->distance(r, u, coincident)};
// Check if this distance is the new minimum.
if (d < min_dist) {
@ -354,7 +354,7 @@ Cell::to_hdf5(hid_t cell_group) const
//==============================================================================
bool
Cell::contains_simple(Position r, Angle a, int32_t on_surface) const
Cell::contains_simple(Position r, Direction u, int32_t on_surface) const
{
for (int32_t token : rpn) {
if (token < OP_UNION) {
@ -367,7 +367,7 @@ Cell::contains_simple(Position r, Angle a, int32_t on_surface) const
return false;
} else {
// Note the off-by-one indexing
bool sense = surfaces_c[abs(token)-1]->sense(r, a);
bool sense = surfaces_c[abs(token)-1]->sense(r, u);
if (sense != (token > 0)) {return false;}
}
}
@ -378,7 +378,7 @@ Cell::contains_simple(Position r, Angle a, int32_t on_surface) const
//==============================================================================
bool
Cell::contains_complex(Position r, Angle a, int32_t on_surface) const
Cell::contains_complex(Position r, Direction u, int32_t on_surface) const
{
// Make a stack of booleans. We don't know how big it needs to be, but we do
// know that rpn.size() is an upper-bound.
@ -409,7 +409,7 @@ Cell::contains_complex(Position r, Angle a, int32_t on_surface) const
stack[i_stack] = false;
} else {
// Note the off-by-one indexing
bool sense = surfaces_c[abs(token)-1]->sense(r, a);;
bool sense = surfaces_c[abs(token)-1]->sense(r, u);;
stack[i_stack] = (sense == (token > 0));
}
}
@ -492,16 +492,16 @@ extern "C" {
bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface)
{
Position r {xyz};
Angle a {uvw};
return c->contains(r, a, on_surface);
Direction u {uvw};
return c->contains(r, u, on_surface);
}
void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface,
double *min_dist, int32_t *i_surf)
{
Position r {xyz};
Angle a {uvw};
std::pair<double, int32_t> out = c->distance(r, a, on_surface);
Direction u {uvw};
std::pair<double, int32_t> out = c->distance(r, u, on_surface);
*min_dist = out.first;
*i_surf = out.second;
}

11
src/cell.h
View file

@ -8,7 +8,8 @@
#include "hdf5.h"
#include "pugixml.hpp"
#include "geometry.h"
#include "position.h"
namespace openmc {
@ -98,19 +99,19 @@ public:
//! known to be on. This index takes precedence over surface sense
//! calculations.
bool
contains(Position r, Angle a, int32_t on_surface) const;
contains(Position r, Direction u, int32_t on_surface) const;
//! Find the oncoming boundary of this cell.
std::pair<double, int32_t>
distance(Position r, Angle a, int32_t on_surface) const;
distance(Position r, Direction u, int32_t on_surface) const;
//! \brief Write cell information to an HDF5 group.
//! @param group_id An HDF5 group id.
void to_hdf5(hid_t group_id) const;
protected:
bool contains_simple(Position r, Angle a, int32_t on_surface) const;
bool contains_complex(Position r, Angle a, int32_t on_surface) const;
bool contains_simple(Position r, Direction u, int32_t on_surface) const;
bool contains_complex(Position r, Direction u, int32_t on_surface) const;
};
} // namespace openmc

64
src/geometry.h
View file

@ -1,68 +1,10 @@
#ifndef GEOMETRY_H
#define GEOMETRY_H
#ifndef OPENMC_GEOMETRY_H
#define OPENMC_GEOMETRY_H
namespace openmc {
extern "C" int openmc_root_universe;
struct Position {
double x = 0.;
double y = 0.;
double z = 0.;
Position() = default;
Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { };
Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
Position& operator+=(Position);
Position& operator+=(double);
Position& operator-=(Position);
Position& operator-=(double);
Position& operator*=(Position);
Position& operator*=(double);
const double& operator[](int i) const {
switch (i) {
case 0: return x;
case 1: return y;
case 2: return z;
}
}
double& operator[](int i) {
switch (i) {
case 0: return x;
case 1: return y;
case 2: return z;
}
}
inline double dot(Position other) {
return x*other.x + y*other.y + z*other.z;
}
};
inline Position operator+(Position a, Position b) { return a += b; }
inline Position operator+(Position a, double b) { return a += b; }
inline Position operator+(double a, Position b) { return b += a; }
inline Position operator-(Position a, Position b) { return a -= b; }
inline Position operator-(Position a, double b) { return a -= b; }
inline Position operator-(double a, Position b) { return b -= a; }
inline Position operator*(Position a, Position b) { return a *= b; }
inline Position operator*(Position a, double b) { return a *= b; }
inline Position operator*(double a, Position b) { return b *= a; }
struct Angle : Position {
double& u() { return x; }
double& v() { return y; }
double& w() { return z; }
Angle() = default;
Angle(double u, double v, double w) : Position{u, v, w} { };
Angle(const double uvw[]) : Position{uvw} { };
Angle(Position r) : Position{r} { };
};
} // namespace openmc
#endif // GEOMETRY_H
#endif // OPENMC_GEOMETRY_H

2
src/hdf5_interface.h
View file

@ -11,7 +11,7 @@
#include <vector>
#include <complex.h>
#include "geometry.h"
#include "position.h"
namespace openmc {

58
src/lattice.cpp
View file

@ -223,25 +223,23 @@ RectLattice::are_valid_indices(const int i_xyz[3]) const
//==============================================================================
std::pair<double, std::array<int, 3>>
RectLattice::distance(Position r, Angle a, const int i_xyz[3]) const
RectLattice::distance(Position r, Direction u, const int i_xyz[3]) const
{
// Get short aliases to the coordinates.
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)};
double y0 {copysign(0.5 * pitch[1], v)};
double x0 {copysign(0.5 * pitch[0], u.x)};
double y0 {copysign(0.5 * pitch[1], u.y)};
// 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) {
if ((std::abs(x - x0) > FP_PRECISION) && u.x != 0) {
d = (x0 - x) / u.x;
if (u.x > 0) {
lattice_trans = {1, 0, 0};
} else {
lattice_trans = {-1, 0, 0};
@ -249,11 +247,11 @@ RectLattice::distance(Position r, Angle a, const int i_xyz[3]) const
}
// Front and back sides
if ((std::abs(y - y0) > FP_PRECISION) && v != 0) {
double this_d = (y0 - y) / v;
if ((std::abs(y - y0) > FP_PRECISION) && u.y != 0) {
double this_d = (y0 - y) / u.y;
if (this_d < d) {
d = this_d;
if (v > 0) {
if (u.y > 0) {
lattice_trans = {0, 1, 0};
} else {
lattice_trans = {0, -1, 0};
@ -263,13 +261,12 @@ RectLattice::distance(Position r, Angle a, const int i_xyz[3]) const
// Top and bottom sides
if (is_3d) {
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;
double z0 {copysign(0.5 * pitch[2], 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 (w > 0) {
if (u.z > 0) {
lattice_trans = {0, 0, 1};
} else {
lattice_trans = {0, 0, -1};
@ -579,11 +576,11 @@ HexLattice::are_valid_indices(const int i_xyz[3]) const
//==============================================================================
std::pair<double, std::array<int, 3>>
HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const
HexLattice::distance(Position r, Direction u, const int i_xyz[3]) const
{
// Compute the direction on the hexagonal basis.
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;
double beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0;
double gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0;
// Note that hexagonal lattice distance calculations are performed
// using the particle's coordinates relative to the neighbor lattice
@ -636,18 +633,18 @@ HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const
}
// Upper and lower sides.
edge = -copysign(0.5*pitch[0], a.v());
if (a.v() > 0) {
edge = -copysign(0.5*pitch[0], u.y);
if (u.y > 0) {
const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]};
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]};
r_t = get_local_position(r, i_xyz_t);
}
if ((std::abs(r_t.y - edge) > FP_PRECISION) && a.v() != 0) {
double this_d = (edge - r_t.y) / a.v();
if ((std::abs(r_t.y - edge) > FP_PRECISION) && u.y != 0) {
double this_d = (edge - r_t.y) / u.y;
if (this_d < d) {
if (a.v() > 0) {
if (u.y > 0) {
lattice_trans = {0, 1, 0};
} else {
lattice_trans = {0, -1, 0};
@ -659,13 +656,12 @@ HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const
// Top and bottom sides
if (is_3d) {
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;
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 (w > 0) {
if (u.z > 0) {
lattice_trans = {0, 0, 1};
} else {
lattice_trans = {0, 0, -1};
@ -897,8 +893,8 @@ extern "C" {
const int i_xyz[3], double *d, int lattice_trans[3])
{
Position r {xyz};
Angle a {uvw};
std::pair<double, std::array<int, 3>> ld {lat->distance(r, a, i_xyz)};
Direction u {uvw};
std::pair<double, std::array<int, 3>> ld {lat->distance(r, u, i_xyz)};
*d = ld.first;
lattice_trans[0] = ld.second[0];
lattice_trans[1] = ld.second[1];

11
src/lattice.h
View file

@ -7,11 +7,12 @@
#include <unordered_map>
#include <vector>
#include "constants.h"
#include "geometry.h"
#include "hdf5.h"
#include "pugixml.hpp"
#include "constants.h"
#include "position.h"
namespace openmc {
@ -82,7 +83,7 @@ public:
//! @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(Position r, Angle a, const int i_xyz[3]) const
distance(Position r, Direction u, const int i_xyz[3]) const
= 0;
//! \brief Find the lattice tile indices for a given point.
@ -194,7 +195,7 @@ public:
bool are_valid_indices(const int i_xyz[3]) const;
std::pair<double, std::array<int, 3>>
distance(Position r, Angle a, const int i_xyz[3]) const;
distance(Position r, Direction u, const int i_xyz[3]) const;
std::array<int, 3> get_indices(Position r) const;
@ -234,7 +235,7 @@ public:
bool are_valid_indices(const int i_xyz[3]) const;
std::pair<double, std::array<int, 3>>
distance(Position r, Angle a, const int i_xyz[3]) const;
distance(Position r, Direction u, const int i_xyz[3]) const;
std::array<int, 3> get_indices(Position r) const;

6
src/geometry.cpp → src/position.cpp
View file

@ -1,7 +1,11 @@
#include "geometry.h"
#include "position.h"
namespace openmc {
//==============================================================================
// Position implementation
//==============================================================================
Position&
Position::operator+=(Position other)
{

74
src/position.h Normal file
View file

@ -0,0 +1,74 @@
#ifndef OPENMC_POSITION_H
#define OPENMC_POSITION_H
namespace openmc {
//==============================================================================
//! Type representing a position in Cartesian coordinates
//==============================================================================
struct Position {
// Constructors
Position() = default;
Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { };
Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
// Unary operators
Position& operator+=(Position);
Position& operator+=(double);
Position& operator-=(Position);
Position& operator-=(double);
Position& operator*=(Position);
Position& operator*=(double);
const double& operator[](int i) const {
switch (i) {
case 0: return x;
case 1: return y;
case 2: return z;
}
}
double& operator[](int i) {
switch (i) {
case 0: return x;
case 1: return y;
case 2: return z;
}
}
// Other member functions
//! Dot product of two vectors
//! \param[in] other Vector to take dot product with
//! \result Resulting dot product
inline double dot(Position other) {
return x*other.x + y*other.y + z*other.z;
}
// Data members
double x = 0.;
double y = 0.;
double z = 0.;
};
// Binary operators
inline Position operator+(Position a, Position b) { return a += b; }
inline Position operator+(Position a, double b) { return a += b; }
inline Position operator+(double a, Position b) { return b += a; }
inline Position operator-(Position a, Position b) { return a -= b; }
inline Position operator-(Position a, double b) { return a -= b; }
inline Position operator-(double a, Position b) { return b -= a; }
inline Position operator*(Position a, Position b) { return a *= b; }
inline Position operator*(Position a, double b) { return a *= b; }
inline Position operator*(double a, Position b) { return b *= a; }
//==============================================================================
//! Type representing a vector direction in Cartesian coordinates
//==============================================================================
using Direction = Position;
} // namespace openmc
#endif // OPENMC_POSITION_H

184
src/surface.cpp
View file

@ -178,7 +178,7 @@ Surface::Surface(pugi::xml_node surf_node)
}
bool
Surface::sense(Position r, Angle a) const
Surface::sense(Position r, Direction u) const
{
// Evaluate the surface equation at the particle's coordinates to determine
// which side the particle is on.
@ -189,22 +189,22 @@ Surface::sense(Position r, Angle a) const
// Particle may be coincident with this surface. To determine the sense, we
// look at the direction of the particle relative to the surface normal (by
// default in the positive direction) via their dot product.
return a.dot(normal(r)) > 0.0;
return u.dot(normal(r)) > 0.0;
}
return f > 0.0;
}
Angle
Surface::reflect(Position r, Angle a) const
Direction
Surface::reflect(Position r, Direction u) const
{
// Determine projection of direction onto normal and squared magnitude of
// normal.
Angle n = normal(r);
const double projection = n.dot(a);
Direction n = normal(r);
const double projection = n.dot(u);
const double magnitude = n.dot(n);
// Reflect direction according to normal.
return a -= (2.0 * projection / magnitude) * n;
return u -= (2.0 * projection / magnitude) * n;
}
void
@ -257,11 +257,11 @@ PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node)
// The template parameter indicates the axis normal to the plane.
template<int i> double
axis_aligned_plane_distance(Position r, Angle a, bool coincident, double offset)
axis_aligned_plane_distance(Position r, Direction u, bool coincident, double offset)
{
const double f = offset - r[i];
if (coincident or std::abs(f) < FP_COINCIDENT or a[i] == 0.0) return INFTY;
const double d = f / a[i];
if (coincident or std::abs(f) < FP_COINCIDENT or u[i] == 0.0) return INFTY;
const double d = f / u[i];
if (d < 0.0) return INFTY;
return d;
}
@ -281,12 +281,12 @@ double SurfaceXPlane::evaluate(Position r) const
return r.x - x0;
}
double SurfaceXPlane::distance(Position r, Angle a, bool coincident) const
double SurfaceXPlane::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_plane_distance<0>(r, a, coincident, x0);
return axis_aligned_plane_distance<0>(r, u, coincident, x0);
}
Angle SurfaceXPlane::normal(Position r) const
Direction SurfaceXPlane::normal(Position r) const
{
return {1., 0., 0.};
}
@ -299,23 +299,23 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
}
bool SurfaceXPlane::periodic_translate(const PeriodicSurface *other, Position& r,
Angle& a) const
Direction& u) const
{
Angle other_n = other->normal(r);
if (other_n.u() == 1 and other_n.v() == 0 and other_n.w() == 0) {
Direction other_n = other->normal(r);
if (other_n.x == 1 and other_n.y == 0 and other_n.z == 0) {
r.x = x0;
return false;
} else {
// Assume the partner is an YPlane (the only supported partner). Use the
// evaluate function to find y0, then adjust position/angle for rotational
// evaluate function to find y0, then adjust position/Direction for rotational
// symmetry.
double y0 = -other->evaluate({0., 0., 0.});
r.y = r.x - x0 + y0;
r.x = x0;
double u = a.u();
a.u() = -a.v();
a.v() = u;
double ux = u.x;
u.x = -u.y;
u.y = ux;
return true;
}
@ -342,12 +342,12 @@ double SurfaceYPlane::evaluate(Position r) const
return r.y - y0;
}
double SurfaceYPlane::distance(Position r, Angle a, bool coincident) const
double SurfaceYPlane::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_plane_distance<1>(r, a, coincident, y0);
return axis_aligned_plane_distance<1>(r, u, coincident, y0);
}
Angle SurfaceYPlane::normal(Position r) const
Direction SurfaceYPlane::normal(Position r) const
{
return {0., 1., 0.};
}
@ -360,24 +360,24 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
}
bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, Position& r,
Angle& a) const
Direction& u) const
{
Angle other_n = other->normal(r);
if (other_n.u() == 0 and other_n.v() == 1 and other_n.w() == 0) {
Direction other_n = other->normal(r);
if (other_n.x == 0 and other_n.y == 1 and other_n.z == 0) {
// The periodic partner is also aligned along y. Just change the y coord.
r.y = y0;
return false;
} else {
// Assume the partner is an XPlane (the only supported partner). Use the
// evaluate function to find x0, then adjust position/angle for rotational
// evaluate function to find x0, then adjust position/Direction for rotational
// symmetry.
double x0 = -other->evaluate({0., 0., 0.});
r.x = r.y - y0 + x0;
r.y = y0;
double u = a.u();
a.u() = a.v();
a.v() = -u;
double ux = u.x;
u.x = u.y;
u.y = -ux;
return true;
}
@ -404,12 +404,12 @@ double SurfaceZPlane::evaluate(Position r) const
return r.z - z0;
}
double SurfaceZPlane::distance(Position r, Angle a, bool coincident) const
double SurfaceZPlane::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_plane_distance<2>(r, a, coincident, z0);
return axis_aligned_plane_distance<2>(r, u, coincident, z0);
}
Angle SurfaceZPlane::normal(Position r) const
Direction SurfaceZPlane::normal(Position r) const
{
return {0., 0., 1.};
}
@ -422,7 +422,7 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
}
bool SurfaceZPlane::periodic_translate(const PeriodicSurface *other, Position& r,
Angle& a) const
Direction& u) const
{
// Assume the other plane is aligned along z. Just change the z coord.
r.z = z0;
@ -452,10 +452,10 @@ SurfacePlane::evaluate(Position r) const
}
double
SurfacePlane::distance(Position r, Angle a, bool coincident) const
SurfacePlane::distance(Position r, Direction u, bool coincident) const
{
const double f = A*r.x + B*r.y + C*r.z - D;
const double projection = A*a.u() + B*a.v() + C*a.w();
const double projection = A*u.x + B*u.y + C*u.z;
if (coincident or std::abs(f) < FP_COINCIDENT or projection == 0.0) {
return INFTY;
} else {
@ -465,7 +465,7 @@ SurfacePlane::distance(Position r, Angle a, bool coincident) const
}
}
Angle
Direction
SurfacePlane::normal(Position r) const
{
return {A, B, C};
@ -479,7 +479,7 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const
}
bool SurfacePlane::periodic_translate(const PeriodicSurface *other, Position& r,
Angle& a) const
Direction& u) const
{
// This function assumes the other plane shares this plane's normal direction.
@ -520,7 +520,7 @@ axis_aligned_cylinder_evaluate(Position r, double offset1,
// The other two parameters indicate the other two axes. offset1 and offset2
// should correspond with i2 and i3, respectively.
template<int i1, int i2, int i3> double
axis_aligned_cylinder_distance(Position r, Angle u,
axis_aligned_cylinder_distance(Position r, Direction u,
bool coincident, double offset1, double offset2, double radius)
{
const double a = 1.0 - u[i1]*u[i1]; // u^2 + v^2
@ -565,14 +565,14 @@ axis_aligned_cylinder_distance(Position r, Angle u,
// The first template parameter indicates which axis the cylinder is aligned to.
// The other two parameters indicate the other two axes. offset1 and offset2
// should correspond with i2 and i3, respectively.
template<int i1, int i2, int i3> Angle
template<int i1, int i2, int i3> Direction
axis_aligned_cylinder_normal(Position r, double offset1, double offset2)
{
Angle a;
a[i2] = 2.0 * (r[i2] - offset1);
a[i3] = 2.0 * (r[i3] - offset2);
a[i1] = 0.0;
return a;
Direction u;
u[i2] = 2.0 * (r[i2] - offset1);
u[i3] = 2.0 * (r[i3] - offset2);
u[i1] = 0.0;
return u;
}
//==============================================================================
@ -590,13 +590,13 @@ double SurfaceXCylinder::evaluate(Position r) const
return axis_aligned_cylinder_evaluate<1, 2>(r, y0, z0, this->r);
}
double SurfaceXCylinder::distance(Position r, Angle a, bool coincident) const
double SurfaceXCylinder::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cylinder_distance<0, 1, 2>(r, a, coincident, y0, z0,
return axis_aligned_cylinder_distance<0, 1, 2>(r, u, coincident, y0, z0,
this->r);
}
Angle SurfaceXCylinder::normal(Position r) const
Direction SurfaceXCylinder::normal(Position r) const
{
return axis_aligned_cylinder_normal<0, 1, 2>(r, y0, z0);
}
@ -624,13 +624,13 @@ double SurfaceYCylinder::evaluate(Position r) const
return axis_aligned_cylinder_evaluate<0, 2>(r, x0, z0, this->r);
}
double SurfaceYCylinder::distance(Position r, Angle a, bool coincident) const
double SurfaceYCylinder::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cylinder_distance<1, 0, 2>(r, a, coincident, x0, z0,
return axis_aligned_cylinder_distance<1, 0, 2>(r, u, coincident, x0, z0,
this->r);
}
Angle SurfaceYCylinder::normal(Position r) const
Direction SurfaceYCylinder::normal(Position r) const
{
return axis_aligned_cylinder_normal<1, 0, 2>(r, x0, z0);
}
@ -657,13 +657,13 @@ double SurfaceZCylinder::evaluate(Position r) const
return axis_aligned_cylinder_evaluate<0, 1>(r, x0, y0, this->r);
}
double SurfaceZCylinder::distance(Position r, Angle a, bool coincident) const
double SurfaceZCylinder::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cylinder_distance<2, 0, 1>(r, a, coincident, x0, y0,
return axis_aligned_cylinder_distance<2, 0, 1>(r, u, coincident, x0, y0,
this->r);
}
Angle SurfaceZCylinder::normal(Position r) const
Direction SurfaceZCylinder::normal(Position r) const
{
return axis_aligned_cylinder_normal<2, 0, 1>(r, x0, y0);
}
@ -693,12 +693,12 @@ double SurfaceSphere::evaluate(Position r) const
return x*x + y*y + z*z - this->r*this->r;
}
double SurfaceSphere::distance(Position r, Angle a, bool coincident) const
double SurfaceSphere::distance(Position r, Direction u, bool coincident) const
{
const double x = r.x - x0;
const double y = r.y - y0;
const double z = r.z - z0;
const double k = x*a.u() + y*a.v() + z*a.w();
const double k = x*u.x + y*u.y + z*u.z;
const double c = x*x + y*y + z*z - this->r*this->r;
const double quad = k*k - c;
@ -731,7 +731,7 @@ double SurfaceSphere::distance(Position r, Angle a, bool coincident) const
}
}
Angle SurfaceSphere::normal(Position r) const
Direction SurfaceSphere::normal(Position r) const
{
return {2.0*(r.x - x0), 2.0*(r.y - y0), 2.0*(r.z - z0)};
}
@ -764,7 +764,7 @@ axis_aligned_cone_evaluate(Position r, double offset1,
// The other two parameters indicate the other two axes. offset1, offset2,
// and offset3 should correspond with i1, i2, and i3, respectively.
template<int i1, int i2, int i3> double
axis_aligned_cone_distance(Position r, Angle u,
axis_aligned_cone_distance(Position r, Direction u,
bool coincident, double offset1, double offset2, double offset3,
double radius_sq)
{
@ -817,15 +817,15 @@ axis_aligned_cone_distance(Position r, Angle u,
// The first template parameter indicates which axis the cone is aligned to.
// The other two parameters indicate the other two axes. offset1, offset2,
// and offset3 should correspond with i1, i2, and i3, respectively.
template<int i1, int i2, int i3> Angle
template<int i1, int i2, int i3> Direction
axis_aligned_cone_normal(Position r, double offset1, double offset2,
double offset3, double radius_sq)
{
Angle a;
a[i1] = -2.0 * radius_sq * (r[i1] - offset1);
a[i2] = 2.0 * (r[i2] - offset2);
a[i3] = 2.0 * (r[i3] - offset3);
return a;
Direction u;
u[i1] = -2.0 * radius_sq * (r[i1] - offset1);
u[i2] = 2.0 * (r[i2] - offset2);
u[i3] = 2.0 * (r[i3] - offset3);
return u;
}
//==============================================================================
@ -843,13 +843,13 @@ double SurfaceXCone::evaluate(Position r) const
return axis_aligned_cone_evaluate<0, 1, 2>(r, x0, y0, z0, r_sq);
}
double SurfaceXCone::distance(Position r, Angle a, bool coincident) const
double SurfaceXCone::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cone_distance<0, 1, 2>(r, a, coincident, x0, y0, z0,
return axis_aligned_cone_distance<0, 1, 2>(r, u, coincident, x0, y0, z0,
r_sq);
}
Angle SurfaceXCone::normal(Position r) const
Direction SurfaceXCone::normal(Position r) const
{
return axis_aligned_cone_normal<0, 1, 2>(r, x0, y0, z0, r_sq);
}
@ -876,13 +876,13 @@ double SurfaceYCone::evaluate(Position r) const
return axis_aligned_cone_evaluate<1, 0, 2>(r, y0, x0, z0, r_sq);
}
double SurfaceYCone::distance(Position r, Angle a, bool coincident) const
double SurfaceYCone::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cone_distance<1, 0, 2>(r, a, coincident, y0, x0, z0,
return axis_aligned_cone_distance<1, 0, 2>(r, u, coincident, y0, x0, z0,
r_sq);
}
Angle SurfaceYCone::normal(Position r) const
Direction SurfaceYCone::normal(Position r) const
{
return axis_aligned_cone_normal<1, 0, 2>(r, y0, x0, z0, r_sq);
}
@ -909,13 +909,13 @@ double SurfaceZCone::evaluate(Position r) const
return axis_aligned_cone_evaluate<2, 0, 1>(r, z0, x0, y0, r_sq);
}
double SurfaceZCone::distance(Position r, Angle a, bool coincident) const
double SurfaceZCone::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cone_distance<2, 0, 1>(r, a, coincident, z0, x0, y0,
return axis_aligned_cone_distance<2, 0, 1>(r, u, coincident, z0, x0, y0,
r_sq);
}
Angle SurfaceZCone::normal(Position r) const
Direction SurfaceZCone::normal(Position r) const
{
return axis_aligned_cone_normal<2, 0, 1>(r, z0, x0, y0, r_sq);
}
@ -949,14 +949,14 @@ SurfaceQuadric::evaluate(Position r) const
}
double
SurfaceQuadric::distance(Position r, Angle ang, bool coincident) const
SurfaceQuadric::distance(Position r, Direction ang, bool coincident) const
{
const double &x = r.x;
const double &y = r.y;
const double &z = r.z;
const double &u = ang.u();
const double &v = ang.v();
const double &w = ang.w();
const double &u = ang.x;
const double &v = ang.y;
const double &w = ang.z;
const double a = A*u*u + B*v*v + C*w*w + D*u*v + E*v*w + F*u*w;
const double k = (A*u*x + B*v*y + C*w*z + 0.5*(D*(u*y + v*x) +
@ -1002,7 +1002,7 @@ SurfaceQuadric::distance(Position r, Angle ang, bool coincident) const
return d;
}
Angle
Direction
SurfaceQuadric::normal(Position r) const
{
const double &x = r.x;
@ -1222,21 +1222,21 @@ extern "C" {
void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
{
Position r {xyz};
Angle a {uvw};
a = surf->reflect(r, a);
Direction u {uvw};
u = surf->reflect(r, u);
uvw[0] = a.u();
uvw[1] = a.v();
uvw[2] = a.w();
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
}
void surface_normal(Surface *surf, double xyz[3], double uvw[3])
{
Position r {xyz};
Angle a = surf->normal(r);
uvw[0] = a.u();
uvw[1] = a.v();
uvw[2] = a.w();
Direction u = surf->normal(r);
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
}
void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);}
@ -1248,16 +1248,16 @@ extern "C" {
double uvw[3])
{
Position r {xyz};
Angle a {uvw};
bool rotational = surf->periodic_translate(other, r, a);
Direction u {uvw};
bool rotational = surf->periodic_translate(other, r, u);
// Copy back to arrays
xyz[0] = r.x;
xyz[1] = r.y;
xyz[2] = r.z;
uvw[0] = a.u();
uvw[1] = a.v();
uvw[2] = a.w();
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
return rotational;
}

94
src/surface.h
View file

@ -9,7 +9,7 @@
#include "pugixml.hpp"
#include "constants.h"
#include "geometry.h"
#include "position.h"
namespace openmc {
@ -71,14 +71,14 @@ public:
//! point is very close to the surface.
//! @return true if the point is on the "positive" side of the surface and
//! false otherwise.
bool sense(Position r, Angle a) const;
bool sense(Position r, Direction u) const;
//! Determine the direction of a ray reflected from the surface.
//! @param r The point at which the ray is incident.
//! @param o A direction. This is both an input and an output parameter.
//! It specifies the icident direction on input and the reflected direction
//! on output.
Angle reflect(Position r, Angle a) const;
Direction reflect(Position r, Direction u) const;
//! Evaluate the equation describing the surface.
//!
@ -92,12 +92,12 @@ public:
//! @param o The direction of the ray.
//! @param coincident A hint to the code that the given point should lie
//! exactly on the surface.
virtual double distance(Position r, Angle a, bool coincident) const = 0;
virtual double distance(Position r, Direction u, bool coincident) const = 0;
//! Compute the local outward normal direction of the surface.
//! @param r A 3D Cartesian coordinate.
//! @return Normal direction
virtual Angle normal(Position r) const = 0;
virtual Direction normal(Position r) const = 0;
//! Write all information needed to reconstruct the surface to an HDF5 group.
//! @param group_id An HDF5 group id.
@ -132,7 +132,7 @@ public:
//! @return true if this surface and its partner make a rotationally-periodic
//! boundary condition.
virtual bool periodic_translate(const PeriodicSurface *other, Position& r,
Angle& a) const = 0;
Direction& u) const = 0;
//! Get the bounding box for this surface.
virtual BoundingBox bounding_box() const = 0;
@ -150,10 +150,10 @@ class SurfaceXPlane : public PeriodicSurface
public:
explicit SurfaceXPlane(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
BoundingBox bounding_box() const;
};
@ -170,10 +170,10 @@ class SurfaceYPlane : public PeriodicSurface
public:
explicit SurfaceYPlane(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
BoundingBox bounding_box() const;
};
@ -190,10 +190,10 @@ class SurfaceZPlane : public PeriodicSurface
public:
explicit SurfaceZPlane(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
BoundingBox bounding_box() const;
};
@ -210,10 +210,10 @@ class SurfacePlane : public PeriodicSurface
public:
explicit SurfacePlane(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
BoundingBox bounding_box() const;
};
@ -231,8 +231,8 @@ class SurfaceXCylinder : public Surface
public:
explicit SurfaceXCylinder(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -249,8 +249,8 @@ class SurfaceYCylinder : public Surface
public:
explicit SurfaceYCylinder(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -267,8 +267,8 @@ class SurfaceZCylinder : public Surface
public:
explicit SurfaceZCylinder(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -285,8 +285,8 @@ class SurfaceSphere : public Surface
public:
explicit SurfaceSphere(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -303,8 +303,8 @@ class SurfaceXCone : public Surface
public:
explicit SurfaceXCone(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -321,8 +321,8 @@ class SurfaceYCone : public Surface
public:
explicit SurfaceYCone(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -339,8 +339,8 @@ class SurfaceZCone : public Surface
public:
explicit SurfaceZCone(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -357,8 +357,8 @@ class SurfaceQuadric : public Surface
public:
explicit SurfaceQuadric(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Angle a, bool coincident) const;
Angle normal(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -367,19 +367,19 @@ public:
//==============================================================================
extern "C" {
Surface* surface_pointer(int surf_ind);
int surface_id(Surface *surf);
int surface_bc(Surface *surf);
bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
double surface_distance(Surface *surf, double xyz[3], double uvw[3],
bool coincident);
void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
void surface_to_hdf5(Surface *surf, hid_t group);
int surface_i_periodic(PeriodicSurface *surf);
bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
double xyz[3], double uvw[3]);
void free_memory_surfaces_c();
Surface* surface_pointer(int surf_ind);
int surface_id(Surface *surf);
int surface_bc(Surface *surf);
bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
double surface_distance(Surface *surf, double xyz[3], double uvw[3],
bool coincident);
void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
void surface_to_hdf5(Surface *surf, hid_t group);
int surface_i_periodic(PeriodicSurface *surf);
bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
double xyz[3], double uvw[3]);
void free_memory_surfaces_c();
}
} // namespace openmc