From d9a7db16777eb6239bcefd1fccf01161eb3fda05 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Jun 2018 21:00:59 -0500 Subject: [PATCH] Introduce Position/Angle types --- CMakeLists.txt | 1 + src/cell.cpp | 34 ++-- src/cell.h | 11 +- src/geometry.cpp | 59 +++++++ src/geometry.h | 62 +++++++ src/surface.cpp | 449 +++++++++++++++++++++++------------------------ src/surface.h | 130 +++++++------- 7 files changed, 421 insertions(+), 325 deletions(-) create mode 100644 src/geometry.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index ea9bd5969..0be9ccb32 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -386,6 +386,7 @@ 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 diff --git a/src/cell.cpp b/src/cell.cpp index 263191b31..972540fbd 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -279,33 +279,31 @@ Cell::Cell(pugi::xml_node cell_node) //============================================================================== bool -Cell::contains(const double xyz[3], const double uvw[3], - int32_t on_surface) const +Cell::contains(Position r, Angle a, int32_t on_surface) const { if (simple) { - return contains_simple(xyz, uvw, on_surface); + return contains_simple(r, a, on_surface); } else { - return contains_complex(xyz, uvw, on_surface); + return contains_complex(r, a, on_surface); } } //============================================================================== std::pair -Cell::distance(const double xyz[3], const double uvw[3], - int32_t on_surface) const +Cell::distance(Position r, Angle a, int32_t on_surface) const { double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; for (int32_t token : rpn) { // Ignore this token if it corresponds to an operator rather than a region. - if (token >= OP_UNION) {continue;} + if (token >= OP_UNION) continue; // 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(xyz, uvw, coincident)}; + double d {surfaces_c[abs(token)-1]->distance(r, a, coincident)}; // Check if this distance is the new minimum. if (d < min_dist) { @@ -356,8 +354,7 @@ Cell::to_hdf5(hid_t cell_group) const //============================================================================== bool -Cell::contains_simple(const double xyz[3], const double uvw[3], - int32_t on_surface) const +Cell::contains_simple(Position r, Angle a, int32_t on_surface) const { for (int32_t token : rpn) { if (token < OP_UNION) { @@ -370,7 +367,7 @@ Cell::contains_simple(const double xyz[3], const double uvw[3], return false; } else { // Note the off-by-one indexing - bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw); + bool sense = surfaces_c[abs(token)-1]->sense(r, a); if (sense != (token > 0)) {return false;} } } @@ -381,8 +378,7 @@ Cell::contains_simple(const double xyz[3], const double uvw[3], //============================================================================== bool -Cell::contains_complex(const double xyz[3], const double uvw[3], - int32_t on_surface) const +Cell::contains_complex(Position r, Angle a, 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. @@ -413,7 +409,7 @@ Cell::contains_complex(const double xyz[3], const double uvw[3], stack[i_stack] = false; } else { // Note the off-by-one indexing - bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw);; + bool sense = surfaces_c[abs(token)-1]->sense(r, a);; stack[i_stack] = (sense == (token > 0)); } } @@ -494,12 +490,18 @@ extern "C" { bool cell_simple(Cell *c) {return c->simple;} bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface) - {return c->contains(xyz, uvw, on_surface);} + { + Position r {xyz}; + Angle a {uvw}; + return c->contains(r, a, on_surface); + } void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface, double *min_dist, int32_t *i_surf) { - std::pair out = c->distance(xyz, uvw, on_surface); + Position r {xyz}; + Angle a {uvw}; + std::pair out = c->distance(r, a, on_surface); *min_dist = out.first; *i_surf = out.second; } diff --git a/src/cell.h b/src/cell.h index 24413b4fd..08a6abf3f 100644 --- a/src/cell.h +++ b/src/cell.h @@ -8,6 +8,7 @@ #include "hdf5.h" #include "pugixml.hpp" +#include "geometry.h" namespace openmc { @@ -97,21 +98,19 @@ public: //! known to be on. This index takes precedence over surface sense //! calculations. bool - contains(const double xyz[3], const double uvw[3], int32_t on_surface) const; + contains(Position r, Angle a, int32_t on_surface) const; //! Find the oncoming boundary of this cell. std::pair - distance(const double xyz[3], const double uvw[3], int32_t on_surface) const; + distance(Position r, Angle a, 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(const double xyz[3], const double uvw[3], - int32_t on_surface) const; - bool contains_complex(const double xyz[3], const double uvw[3], - int32_t on_surface) const; + bool contains_simple(Position r, Angle a, int32_t on_surface) const; + bool contains_complex(Position r, Angle a, int32_t on_surface) const; }; } // namespace openmc diff --git a/src/geometry.cpp b/src/geometry.cpp new file mode 100644 index 000000000..53bf32b5c --- /dev/null +++ b/src/geometry.cpp @@ -0,0 +1,59 @@ +#include "geometry.h" + +namespace openmc { + +Position& +Position::operator+=(Position other) +{ + x += other.x; + y += other.y; + z += other.z; + return *this; +} + +Position& +Position::operator+=(double v) +{ + x += v; + y += v; + z += v; + return *this; +} + +Position& +Position::operator-=(Position other) +{ + x -= other.x; + y -= other.y; + z -= other.z; + return *this; +} + +Position& +Position::operator-=(double v) +{ + x -= v; + y -= v; + z -= v; + return *this; +} + +Position& +Position::operator*=(Position other) +{ + x *= other.x; + y *= other.y; + z *= other.z; + return *this; +} + +Position& +Position::operator*=(double v) +{ + x *= v; + y *= v; + z *= v; + return *this; +} + +} // namespace openmc diff --git a/src/geometry.h b/src/geometry.h index 8837c961a..792d653ae 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -1,6 +1,68 @@ #ifndef GEOMETRY_H #define 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(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(double uvw[]) : Position{uvw} { }; + Angle(Position r) : Position{r} { }; +}; + +} // namespace openmc + #endif // GEOMETRY_H diff --git a/src/surface.cpp b/src/surface.cpp index 6400da2f7..aec0972b8 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include "error.h" #include "hdf5_interface.h" @@ -177,38 +178,33 @@ Surface::Surface(pugi::xml_node surf_node) } bool -Surface::sense(const double xyz[3], const double uvw[3]) const +Surface::sense(Position r, Angle a) const { // Evaluate the surface equation at the particle's coordinates to determine // which side the particle is on. - const double f = evaluate(xyz); + const double f = evaluate(r); // Check which side of surface the point is on. if (std::abs(f) < FP_COINCIDENT) { // 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. - double norm[3]; - normal(xyz, norm); - return uvw[0] * norm[0] + uvw[1] * norm[1] + uvw[2] * norm[2] > 0.0; + return a.dot(normal(r)) > 0.0; } return f > 0.0; } -void -Surface::reflect(const double xyz[3], double uvw[3]) const +Angle +Surface::reflect(Position r, Angle a) const { // Determine projection of direction onto normal and squared magnitude of // normal. - double norm[3]; - normal(xyz, norm); - const double projection = norm[0]*uvw[0] + norm[1]*uvw[1] + norm[2]*uvw[2]; - const double magnitude = norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2]; + Angle n = normal(r); + const double projection = n.dot(a); + const double magnitude = n.dot(n); // Reflect direction according to normal. - uvw[0] -= 2.0 * projection / magnitude * norm[0]; - uvw[1] -= 2.0 * projection / magnitude * norm[1]; - uvw[2] -= 2.0 * projection / magnitude * norm[2]; + return a -= (2.0 * projection / magnitude) * n; } void @@ -261,33 +257,15 @@ PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node) // The template parameter indicates the axis normal to the plane. template double -axis_aligned_plane_evaluate(const double xyz[3], double offset) +axis_aligned_plane_distance(Position r, Angle a, bool coincident, double offset) { - return xyz[i] - offset; -} - -// The template parameter indicates the axis normal to the plane. -template double -axis_aligned_plane_distance(const double xyz[3], const double uvw[3], - bool coincident, double offset) -{ - const double f = offset - xyz[i]; - if (coincident or std::abs(f) < FP_COINCIDENT or uvw[i] == 0.0) return INFTY; - const double d = f / uvw[i]; + 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 (d < 0.0) return INFTY; return d; } -// The first template parameter indicates the axis normal to the plane. The -// other two parameters indicate the other two axes. -template void -axis_aligned_plane_normal(const double xyz[3], double uvw[3]) -{ - uvw[i1] = 1.0; - uvw[i2] = 0.0; - uvw[i3] = 0.0; -} - //============================================================================== // SurfaceXPlane implementation //============================================================================== @@ -298,20 +276,19 @@ SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0); } -inline double SurfaceXPlane::evaluate(const double xyz[3]) const +double SurfaceXPlane::evaluate(Position r) const { - return axis_aligned_plane_evaluate<0>(xyz, x0); + return r.x - x0; } -inline double SurfaceXPlane::distance(const double xyz[3], const double uvw[3], - bool coincident) const +double SurfaceXPlane::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_plane_distance<0>(xyz, uvw, coincident, x0); + return axis_aligned_plane_distance<0>(r, a, coincident, x0); } -inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceXPlane::normal(Position r) const { - axis_aligned_plane_normal<0, 1, 2>(xyz, uvw); + return {1., 0., 0.}; } void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const @@ -321,26 +298,24 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], - double uvw[3]) const +bool SurfaceXPlane::periodic_translate(const PeriodicSurface *other, Position& r, + Angle& a) const { - double other_norm[3]; - other->normal(xyz, other_norm); - if (other_norm[0] == 1 and other_norm[1] == 0 and other_norm[2] == 0) { - xyz[0] = x0; + Angle other_n = other->normal(r); + if (other_n.u() == 1 and other_n.v() == 0 and other_n.w() == 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 xyz and uvw for rotational + // evaluate function to find y0, then adjust position/angle for rotational // symmetry. - double xyz_test[3] {0, 0, 0}; - double y0 = -other->evaluate(xyz_test); - xyz[1] = xyz[0] - x0 + y0; - xyz[0] = x0; + double y0 = -other->evaluate({0., 0., 0.}); + r.y = r.x - x0 + y0; + r.x = x0; - double u = uvw[0]; - uvw[0] = -uvw[1]; - uvw[1] = u; + double u = a.u(); + a.u() = -a.v(); + a.v() = u; return true; } @@ -349,8 +324,7 @@ bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], BoundingBox SurfaceXPlane::bounding_box() const { - BoundingBox out {x0, x0, -INFTY, INFTY, -INFTY, INFTY}; - return out; + return {x0, x0, -INFTY, INFTY, -INFTY, INFTY}; } //============================================================================== @@ -363,20 +337,19 @@ SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) read_coeffs(surf_node, id, y0); } -inline double SurfaceYPlane::evaluate(const double xyz[3]) const +double SurfaceYPlane::evaluate(Position r) const { - return axis_aligned_plane_evaluate<1>(xyz, y0); + return r.y - y0; } -inline double SurfaceYPlane::distance(const double xyz[3], const double uvw[3], - bool coincident) const +double SurfaceYPlane::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_plane_distance<1>(xyz, uvw, coincident, y0); + return axis_aligned_plane_distance<1>(r, a, coincident, y0); } -inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceYPlane::normal(Position r) const { - axis_aligned_plane_normal<1, 0, 2>(xyz, uvw); + return {0., 1., 0.}; } void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const @@ -386,27 +359,25 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], - double uvw[3]) const +bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, Position& r, + Angle& a) const { - double other_norm[3]; - other->normal(xyz, other_norm); - if (other_norm[0] == 0 and other_norm[1] == 1 and other_norm[2] == 0) { + Angle other_n = other->normal(r); + if (other_n.u() == 0 and other_n.v() == 1 and other_n.w() == 0) { // The periodic partner is also aligned along y. Just change the y coord. - xyz[1] = y0; + 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 xyz and uvw for rotational + // evaluate function to find x0, then adjust position/angle for rotational // symmetry. - double xyz_test[3] {0, 0, 0}; - double x0 = -other->evaluate(xyz_test); - xyz[0] = xyz[1] - y0 + x0; - xyz[1] = y0; + double x0 = -other->evaluate({0., 0., 0.}); + r.x = r.y - y0 + x0; + r.y = y0; - double u = uvw[0]; - uvw[0] = uvw[1]; - uvw[1] = -u; + double u = a.u(); + a.u() = a.v(); + a.v() = -u; return true; } @@ -415,8 +386,7 @@ bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], BoundingBox SurfaceYPlane::bounding_box() const { - BoundingBox out {-INFTY, INFTY, y0, y0, -INFTY, INFTY}; - return out; + return {-INFTY, INFTY, y0, y0, -INFTY, INFTY}; } //============================================================================== @@ -429,20 +399,19 @@ SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) read_coeffs(surf_node, id, z0); } -inline double SurfaceZPlane::evaluate(const double xyz[3]) const +double SurfaceZPlane::evaluate(Position r) const { - return axis_aligned_plane_evaluate<2>(xyz, z0); + return r.z - z0; } -inline double SurfaceZPlane::distance(const double xyz[3], const double uvw[3], - bool coincident) const +double SurfaceZPlane::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_plane_distance<2>(xyz, uvw, coincident, z0); + return axis_aligned_plane_distance<2>(r, a, coincident, z0); } -inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceZPlane::normal(Position r) const { - axis_aligned_plane_normal<2, 0, 1>(xyz, uvw); + return {0., 0., 1.}; } void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const @@ -452,19 +421,18 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3], - double uvw[3]) const +bool SurfaceZPlane::periodic_translate(const PeriodicSurface *other, Position& r, + Angle& a) const { // Assume the other plane is aligned along z. Just change the z coord. - xyz[2] = z0; + r.z = z0; return false; } BoundingBox SurfaceZPlane::bounding_box() const { - BoundingBox out {-INFTY, INFTY, -INFTY, INFTY, z0, z0}; - return out; + return {-INFTY, INFTY, -INFTY, INFTY, z0, z0}; } //============================================================================== @@ -478,17 +446,16 @@ SurfacePlane::SurfacePlane(pugi::xml_node surf_node) } double -SurfacePlane::evaluate(const double xyz[3]) const +SurfacePlane::evaluate(Position r) const { - return A*xyz[0] + B*xyz[1] + C*xyz[2] - D; + return A*r.x + B*r.y + C*r.z - D; } double -SurfacePlane::distance(const double xyz[3], const double uvw[3], - bool coincident) const +SurfacePlane::distance(Position r, Angle a, bool coincident) const { - const double f = A*xyz[0] + B*xyz[1] + C*xyz[2] - D; - const double projection = A*uvw[0] + B*uvw[1] + C*uvw[2]; + 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(); if (coincident or std::abs(f) < FP_COINCIDENT or projection == 0.0) { return INFTY; } else { @@ -498,12 +465,10 @@ SurfacePlane::distance(const double xyz[3], const double uvw[3], } } -void -SurfacePlane::normal(const double xyz[3], double uvw[3]) const +Angle +SurfacePlane::normal(Position r) const { - uvw[0] = A; - uvw[1] = B; - uvw[2] = C; + return {A, B, C}; } void SurfacePlane::to_hdf5_inner(hid_t group_id) const @@ -513,18 +478,18 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3], - double uvw[3]) const +bool SurfacePlane::periodic_translate(const PeriodicSurface *other, Position& r, + Angle& a) const { // This function assumes the other plane shares this plane's normal direction. // Determine the distance to intersection. - double d = evaluate(xyz) / (A*A + B*B + C*C); + double d = evaluate(r) / (A*A + B*B + C*C); // Move the particle that distance along the normal vector. - xyz[0] -= d * A; - xyz[1] -= d * B; - xyz[2] -= d * C; + r.x -= d * A; + r.y -= d * B; + r.z -= d * C; return false; } @@ -532,8 +497,7 @@ bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3], BoundingBox SurfacePlane::bounding_box() const { - BoundingBox out {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; - return out; + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; } //============================================================================== @@ -544,28 +508,28 @@ SurfacePlane::bounding_box() const // cylinder. offset1 and offset2 should correspond with i1 and i2, // respectively. template double -axis_aligned_cylinder_evaluate(const double xyz[3], double offset1, +axis_aligned_cylinder_evaluate(Position r, double offset1, double offset2, double radius) { - const double xyz1 = xyz[i1] - offset1; - const double xyz2 = xyz[i2] - offset2; - return xyz1*xyz1 + xyz2*xyz2 - radius*radius; + const double r1 = r[i1] - offset1; + const double r2 = r[i2] - offset2; + return r1*r1 + r2*r2 - radius*radius; } // 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 double -axis_aligned_cylinder_distance(const double xyz[3], const double uvw[3], +axis_aligned_cylinder_distance(Position r, Angle u, bool coincident, double offset1, double offset2, double radius) { - const double a = 1.0 - uvw[i1]*uvw[i1]; // u^2 + v^2 + const double a = 1.0 - u[i1]*u[i1]; // u^2 + v^2 if (a == 0.0) return INFTY; - const double xyz2 = xyz[i2] - offset1; - const double xyz3 = xyz[i3] - offset2; - const double k = xyz2 * uvw[i2] + xyz3 * uvw[i3]; - const double c = xyz2*xyz2 + xyz3*xyz3 - radius*radius; + const double r2 = r[i2] - offset1; + const double r3 = r[i3] - offset2; + const double k = r2 * u[i2] + r3 * u[i3]; + const double c = r2*r2 + r3*r3 - radius*radius; const double quad = k*k - a*c; if (quad < 0.0) { @@ -601,13 +565,14 @@ axis_aligned_cylinder_distance(const double xyz[3], const double uvw[3], // 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 void -axis_aligned_cylinder_normal(const double xyz[3], double uvw[3], double offset1, - double offset2) +template Angle +axis_aligned_cylinder_normal(Position r, double offset1, double offset2) { - uvw[i2] = 2.0 * (xyz[i2] - offset1); - uvw[i3] = 2.0 * (xyz[i3] - offset2); - uvw[i1] = 0.0; + Angle a; + a[i2] = 2.0 * (r[i2] - offset1); + a[i3] = 2.0 * (r[i3] - offset2); + a[i1] = 0.0; + return a; } //============================================================================== @@ -620,21 +585,20 @@ SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node) read_coeffs(surf_node, id, y0, z0, r); } -inline double SurfaceXCylinder::evaluate(const double xyz[3]) const +double SurfaceXCylinder::evaluate(Position r) const { - return axis_aligned_cylinder_evaluate<1, 2>(xyz, y0, z0, r); + return axis_aligned_cylinder_evaluate<1, 2>(r, y0, z0, this->r); } -inline double SurfaceXCylinder::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceXCylinder::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cylinder_distance<0, 1, 2>(xyz, uvw, coincident, y0, z0, - r); + return axis_aligned_cylinder_distance<0, 1, 2>(r, a, coincident, y0, z0, + this->r); } -inline void SurfaceXCylinder::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceXCylinder::normal(Position r) const { - axis_aligned_cylinder_normal<0, 1, 2>(xyz, uvw, y0, z0); + return axis_aligned_cylinder_normal<0, 1, 2>(r, y0, z0); } @@ -655,21 +619,20 @@ SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, z0, r); } -inline double SurfaceYCylinder::evaluate(const double xyz[3]) const +double SurfaceYCylinder::evaluate(Position r) const { - return axis_aligned_cylinder_evaluate<0, 2>(xyz, x0, z0, r); + return axis_aligned_cylinder_evaluate<0, 2>(r, x0, z0, this->r); } -inline double SurfaceYCylinder::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceYCylinder::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cylinder_distance<1, 0, 2>(xyz, uvw, coincident, x0, z0, - r); + return axis_aligned_cylinder_distance<1, 0, 2>(r, a, coincident, x0, z0, + this->r); } -inline void SurfaceYCylinder::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceYCylinder::normal(Position r) const { - axis_aligned_cylinder_normal<1, 0, 2>(xyz, uvw, x0, z0); + return axis_aligned_cylinder_normal<1, 0, 2>(r, x0, z0); } void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const @@ -689,21 +652,20 @@ SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, y0, r); } -inline double SurfaceZCylinder::evaluate(const double xyz[3]) const +double SurfaceZCylinder::evaluate(Position r) const { - return axis_aligned_cylinder_evaluate<0, 1>(xyz, x0, y0, r); + return axis_aligned_cylinder_evaluate<0, 1>(r, x0, y0, this->r); } -inline double SurfaceZCylinder::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceZCylinder::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cylinder_distance<2, 0, 1>(xyz, uvw, coincident, x0, y0, - r); + return axis_aligned_cylinder_distance<2, 0, 1>(r, a, coincident, x0, y0, + this->r); } -inline void SurfaceZCylinder::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceZCylinder::normal(Position r) const { - axis_aligned_cylinder_normal<2, 0, 1>(xyz, uvw, x0, y0); + return axis_aligned_cylinder_normal<2, 0, 1>(r, x0, y0); } void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const @@ -723,22 +685,21 @@ SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, y0, z0, r); } -double SurfaceSphere::evaluate(const double xyz[3]) const +double SurfaceSphere::evaluate(Position r) const { - const double x = xyz[0] - x0; - const double y = xyz[1] - y0; - const double z = xyz[2] - z0; - return x*x + y*y + z*z - r*r; + const double x = r.x - x0; + const double y = r.y - y0; + const double z = r.z - z0; + return x*x + y*y + z*z - this->r*this->r; } -double SurfaceSphere::distance(const double xyz[3], const double uvw[3], - bool coincident) const +double SurfaceSphere::distance(Position r, Angle a, bool coincident) const { - const double x = xyz[0] - x0; - const double y = xyz[1] - y0; - const double z = xyz[2] - z0; - const double k = x*uvw[0] + y*uvw[1] + z*uvw[2]; - const double c = x*x + y*y + z*z - r*r; + 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 c = x*x + y*y + z*z - this->r*this->r; const double quad = k*k - c; if (quad < 0.0) { @@ -770,11 +731,9 @@ double SurfaceSphere::distance(const double xyz[3], const double uvw[3], } } -inline void SurfaceSphere::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceSphere::normal(Position r) const { - uvw[0] = 2.0 * (xyz[0] - x0); - uvw[1] = 2.0 * (xyz[1] - y0); - uvw[2] = 2.0 * (xyz[2] - z0); + return {2.0*(r.x - x0), 2.0*(r.y - y0), 2.0*(r.z - z0)}; } void SurfaceSphere::to_hdf5_inner(hid_t group_id) const @@ -792,30 +751,30 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const // The other two parameters indicate the other two axes. offset1, offset2, // and offset3 should correspond with i1, i2, and i3, respectively. template double -axis_aligned_cone_evaluate(const double xyz[3], double offset1, +axis_aligned_cone_evaluate(Position r, double offset1, double offset2, double offset3, double radius_sq) { - const double xyz1 = xyz[i1] - offset1; - const double xyz2 = xyz[i2] - offset2; - const double xyz3 = xyz[i3] - offset3; - return xyz2*xyz2 + xyz3*xyz3 - radius_sq*xyz1*xyz1; + const double r1 = r[i1] - offset1; + const double r2 = r[i2] - offset2; + const double r3 = r[i3] - offset3; + return r2*r2 + r3*r3 - radius_sq*r1*r1; } // 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 double -axis_aligned_cone_distance(const double xyz[3], const double uvw[3], +axis_aligned_cone_distance(Position r, Angle u, bool coincident, double offset1, double offset2, double offset3, double radius_sq) { - const double xyz1 = xyz[i1] - offset1; - const double xyz2 = xyz[i2] - offset2; - const double xyz3 = xyz[i3] - offset3; - const double a = uvw[i2]*uvw[i2] + uvw[i3]*uvw[i3] - - radius_sq*uvw[i1]*uvw[i1]; - const double k = xyz2*uvw[i2] + xyz3*uvw[i3] - radius_sq*xyz1*uvw[i1]; - const double c = xyz2*xyz2 + xyz3*xyz3 - radius_sq*xyz1*xyz1; + const double r1 = r[i1] - offset1; + const double r2 = r[i2] - offset2; + const double r3 = r[i3] - offset3; + const double a = u[i2]*u[i2] + u[i3]*u[i3] + - radius_sq*u[i1]*u[i1]; + const double k = r2*u[i2] + r3*u[i3] - radius_sq*r1*u[i1]; + const double c = r2*r2 + r3*r3 - radius_sq*r1*r1; double quad = k*k - a*c; double d; @@ -858,13 +817,15 @@ axis_aligned_cone_distance(const double xyz[3], const double uvw[3], // 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 void -axis_aligned_cone_normal(const double xyz[3], double uvw[3], double offset1, - double offset2, double offset3, double radius_sq) +template Angle +axis_aligned_cone_normal(Position r, double offset1, double offset2, + double offset3, double radius_sq) { - uvw[i1] = -2.0 * radius_sq * (xyz[i1] - offset1); - uvw[i2] = 2.0 * (xyz[i2] - offset2); - uvw[i3] = 2.0 * (xyz[i3] - offset3); + 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; } //============================================================================== @@ -877,21 +838,20 @@ SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, y0, z0, r_sq); } -inline double SurfaceXCone::evaluate(const double xyz[3]) const +double SurfaceXCone::evaluate(Position r) const { - return axis_aligned_cone_evaluate<0, 1, 2>(xyz, x0, y0, z0, r_sq); + return axis_aligned_cone_evaluate<0, 1, 2>(r, x0, y0, z0, r_sq); } -inline double SurfaceXCone::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceXCone::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cone_distance<0, 1, 2>(xyz, uvw, coincident, x0, y0, z0, + return axis_aligned_cone_distance<0, 1, 2>(r, a, coincident, x0, y0, z0, r_sq); } -inline void SurfaceXCone::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceXCone::normal(Position r) const { - axis_aligned_cone_normal<0, 1, 2>(xyz, uvw, x0, y0, z0, r_sq); + return axis_aligned_cone_normal<0, 1, 2>(r, x0, y0, z0, r_sq); } void SurfaceXCone::to_hdf5_inner(hid_t group_id) const @@ -911,21 +871,20 @@ SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, y0, z0, r_sq); } -inline double SurfaceYCone::evaluate(const double xyz[3]) const +double SurfaceYCone::evaluate(Position r) const { - return axis_aligned_cone_evaluate<1, 0, 2>(xyz, y0, x0, z0, r_sq); + return axis_aligned_cone_evaluate<1, 0, 2>(r, y0, x0, z0, r_sq); } -inline double SurfaceYCone::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceYCone::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cone_distance<1, 0, 2>(xyz, uvw, coincident, y0, x0, z0, + return axis_aligned_cone_distance<1, 0, 2>(r, a, coincident, y0, x0, z0, r_sq); } -inline void SurfaceYCone::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceYCone::normal(Position r) const { - axis_aligned_cone_normal<1, 0, 2>(xyz, uvw, y0, x0, z0, r_sq); + return axis_aligned_cone_normal<1, 0, 2>(r, y0, x0, z0, r_sq); } void SurfaceYCone::to_hdf5_inner(hid_t group_id) const @@ -945,21 +904,20 @@ SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, y0, z0, r_sq); } -inline double SurfaceZCone::evaluate(const double xyz[3]) const +double SurfaceZCone::evaluate(Position r) const { - return axis_aligned_cone_evaluate<2, 0, 1>(xyz, z0, x0, y0, r_sq); + return axis_aligned_cone_evaluate<2, 0, 1>(r, z0, x0, y0, r_sq); } -inline double SurfaceZCone::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceZCone::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cone_distance<2, 0, 1>(xyz, uvw, coincident, z0, x0, y0, + return axis_aligned_cone_distance<2, 0, 1>(r, a, coincident, z0, x0, y0, r_sq); } -inline void SurfaceZCone::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceZCone::normal(Position r) const { - axis_aligned_cone_normal<2, 0, 1>(xyz, uvw, z0, x0, y0, r_sq); + return axis_aligned_cone_normal<2, 0, 1>(r, z0, x0, y0, r_sq); } void SurfaceZCone::to_hdf5_inner(hid_t group_id) const @@ -980,26 +938,25 @@ SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node) } double -SurfaceQuadric::evaluate(const double xyz[3]) const +SurfaceQuadric::evaluate(Position r) const { - const double &x = xyz[0]; - const double &y = xyz[1]; - const double &z = xyz[2]; + const double x = r.x; + const double y = r.y; + const double z = r.z; return x*(A*x + D*y + G) + y*(B*y + E*z + H) + z*(C*z + F*x + J) + K; } double -SurfaceQuadric::distance(const double xyz[3], - const double uvw[3], bool coincident) const +SurfaceQuadric::distance(Position r, Angle ang, bool coincident) const { - const double &x = xyz[0]; - const double &y = xyz[1]; - const double &z = xyz[2]; - const double &u = uvw[0]; - const double &v = uvw[1]; - const double &w = uvw[2]; + 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 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) + @@ -1045,15 +1002,15 @@ SurfaceQuadric::distance(const double xyz[3], return d; } -void -SurfaceQuadric::normal(const double xyz[3], double uvw[3]) const +Angle +SurfaceQuadric::normal(Position r) const { - const double &x = xyz[0]; - const double &y = xyz[1]; - const double &z = xyz[2]; - uvw[0] = 2.0*A*x + D*y + F*z + G; - uvw[1] = 2.0*B*y + D*x + E*z + H; - uvw[2] = 2.0*C*z + E*y + F*x + J; + const double &x = r.x; + const double &y = r.y; + const double &z = r.z; + return {2.0*A*x + D*y + F*z + G, + 2.0*B*y + D*x + E*z + H, + 2.0*C*z + E*y + F*x + J}; } void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const @@ -1263,10 +1220,24 @@ extern "C" { int surface_bc(Surface *surf) {return surf->bc;} void surface_reflect(Surface *surf, double xyz[3], double uvw[3]) - {surf->reflect(xyz, uvw);} + { + Position r {xyz}; + Angle a {uvw}; + a = surf->reflect(r, a); + + uvw[0] = a.u(); + uvw[1] = a.v(); + uvw[2] = a.w(); + } void surface_normal(Surface *surf, double xyz[3], double uvw[3]) - {return surf->normal(xyz, uvw);} + { + Position r {xyz}; + Angle a = surf->normal(r); + uvw[0] = a.u(); + uvw[1] = a.v(); + uvw[2] = a.w(); + } void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);} @@ -1275,7 +1246,21 @@ extern "C" { bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3], double uvw[3]) - {return surf->periodic_translate(other, xyz, uvw);} + { + Position r {xyz}; + Angle a {uvw}; + bool rotational = surf->periodic_translate(other, r, a); + + // 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(); + + return rotational; + } void free_memory_surfaces_c() { diff --git a/src/surface.h b/src/surface.h index 666fc0dba..be3ad52ab 100644 --- a/src/surface.h +++ b/src/surface.h @@ -9,6 +9,7 @@ #include "pugixml.hpp" #include "constants.h" +#include "geometry.h" namespace openmc { @@ -65,39 +66,38 @@ public: virtual ~Surface() {} //! Determine which side of a surface a point lies on. - //! @param xyz[3] The 3D Cartesian coordinate of a point. - //! @param uvw[3] A direction used to "break ties" and pick a sense when the + //! @param r The 3D Cartesian coordinate of a point. + //! @param o A direction used to "break ties" and pick a sense when the //! 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(const double xyz[3], const double uvw[3]) const; + bool sense(Position r, Angle a) const; //! Determine the direction of a ray reflected from the surface. - //! @param xyz[3] The point at which the ray is incident. - //! @param uvw[3] A direction. This is both an input and an output parameter. + //! @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. - void reflect(const double xyz[3], double uvw[3]) const; + Angle reflect(Position r, Angle a) const; //! Evaluate the equation describing the surface. //! //! Surfaces can be described by some function f(x, y, z) = 0. This member //! function evaluates that mathematical function. - //! @param xyz[3] A 3D Cartesian coordinate. - virtual double evaluate(const double xyz[3]) const = 0; + //! @param r A 3D Cartesian coordinate. + virtual double evaluate(Position r) const = 0; //! Compute the distance between a point and the surface along a ray. - //! @param xyz[3] A 3D Cartesian coordinate. - //! @param uvw[3] The direction of the ray. + //! @param r A 3D Cartesian coordinate. + //! @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(const double xyz[3], const double uvw[3], - bool coincident) const = 0; + virtual double distance(Position r, Angle a, bool coincident) const = 0; //! Compute the local outward normal direction of the surface. - //! @param xyz[3] A 3D Cartesian coordinate. - //! @param uvw[3] This output argument provides the normal. - virtual void normal(const double xyz[3], double uvw[3]) const = 0; + //! @param r A 3D Cartesian coordinate. + //! @return Normal direction + virtual Angle normal(Position r) const = 0; //! Write all information needed to reconstruct the surface to an HDF5 group. //! @param group_id An HDF5 group id. @@ -125,14 +125,14 @@ public: //! Translate a particle onto this surface from a periodic partner surface. //! @param other A pointer to the partner surface in this periodic BC. - //! @param xyz[3] A point on the partner surface that will be translated onto + //! @param r A point on the partner surface that will be translated onto //! this surface. - //! @param uvw[3] A direction that will be rotated for systems with rotational + //! @param a A direction that will be rotated for systems with rotational //! periodicity. //! @return true if this surface and its partner make a rotationally-periodic //! boundary condition. - virtual bool periodic_translate(PeriodicSurface *other, double xyz[3], - double uvw[3]) const = 0; + virtual bool periodic_translate(const PeriodicSurface *other, Position& r, + Angle& a) const = 0; //! Get the bounding box for this surface. virtual BoundingBox bounding_box() const = 0; @@ -149,12 +149,11 @@ class SurfaceXPlane : public PeriodicSurface double x0; public: explicit SurfaceXPlane(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], bool coincident) - const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) const; BoundingBox bounding_box() const; }; @@ -170,12 +169,11 @@ class SurfaceYPlane : public PeriodicSurface double y0; public: explicit SurfaceYPlane(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) const; BoundingBox bounding_box() const; }; @@ -191,12 +189,11 @@ class SurfaceZPlane : public PeriodicSurface double z0; public: explicit SurfaceZPlane(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) const; BoundingBox bounding_box() const; }; @@ -212,12 +209,11 @@ class SurfacePlane : public PeriodicSurface double A, B, C, D; public: explicit SurfacePlane(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], bool coincident) - const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) const; BoundingBox bounding_box() const; }; @@ -234,10 +230,9 @@ class SurfaceXCylinder : public Surface double y0, z0, r; public: explicit SurfaceXCylinder(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -253,10 +248,9 @@ class SurfaceYCylinder : public Surface double x0, z0, r; public: explicit SurfaceYCylinder(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -272,10 +266,9 @@ class SurfaceZCylinder : public Surface double x0, y0, r; public: explicit SurfaceZCylinder(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -291,10 +284,9 @@ class SurfaceSphere : public Surface double x0, y0, z0, r; public: explicit SurfaceSphere(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -310,10 +302,9 @@ class SurfaceXCone : public Surface double x0, y0, z0, r_sq; public: explicit SurfaceXCone(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -329,10 +320,9 @@ class SurfaceYCone : public Surface double x0, y0, z0, r_sq; public: explicit SurfaceYCone(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -348,10 +338,9 @@ class SurfaceZCone : public Surface double x0, y0, z0, r_sq; public: explicit SurfaceZCone(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -367,10 +356,9 @@ class SurfaceQuadric : public Surface double A, B, C, D, E, F, G, H, J, K; public: explicit SurfaceQuadric(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; };