From 4cd335f6a93b21d5010e9f9981093a05a5d8435e Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 11 Aug 2017 22:05:33 -0400 Subject: [PATCH] Add x-, y- planes and cylinders to C++ Introudce templated generic functions to reduce code repetition for surfaces --- src/surface_header.C | 358 +++++++++++++++++++++++++++++++++---------- 1 file changed, 276 insertions(+), 82 deletions(-) diff --git a/src/surface_header.C b/src/surface_header.C index 34ac7cf336..f62be77ab3 100644 --- a/src/surface_header.C +++ b/src/surface_header.C @@ -33,36 +33,34 @@ class Surface { char name[104]; // User-defined name public: - bool sense(double xyz[3], double uvw[3]); - void reflect(double xyz[3], double uvw[3]); - virtual double evaluate(double xyz[3]) = 0; - virtual double distance(double xyz[3], double uvw[3], bool coincident) = 0; - virtual void normal(double xyz[3], double uvw[3]) = 0; + bool sense(const double xyz[3], const double uvw[3]) const; + void reflect(const double xyz[3], double uvw[3]) const; + virtual double evaluate(const double xyz[3]) const = 0; + virtual double distance(const double xyz[3], const double uvw[3], + bool coincident) const = 0; + virtual void normal(const double xyz[3], double uvw[3]) const = 0; }; bool -Surface::sense(double xyz[3], double uvw[3]) { +Surface::sense(const double xyz[3], const double uvw[3]) const { // Evaluate the surface equation at the particle's coordinates to determine // which side the particle is on. const double f = evaluate(xyz); // Check which side of surface the point is on. - bool s; if (fabs(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); - s = (uvw[0] * norm[0] + uvw[1] * norm[1] + uvw[2] * norm[2] > 0.0); - } else { - s = (f > 0.0); + return uvw[0] * norm[0] + uvw[1] * norm[1] + uvw[2] * norm[2] > 0.0; } - return s; + return f > 0.0; } void -Surface::reflect(double xyz[3], double uvw[3]) { +Surface::reflect(const double xyz[3], double uvw[3]) const { // Determine projection of direction onto normal and squared magnitude of // normal. double norm[3]; @@ -76,15 +74,113 @@ Surface::reflect(double xyz[3], double uvw[3]) { uvw[2] -= 2.0 * projection / magnitude * norm[2]; } +//============================================================================== +// Generic functions for X-, Y-, and Z-, planes +//============================================================================== + +template double +axis_aligned_plane_evaluate(const double xyz[3], double offset) { + return xyz[i] - offset;} + +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 fabs(f) < FP_COINCIDENT or uvw[i] == 0.0) return INFTY; + const double d = f / uvw[i]; + if (d < 0.0) return INFTY; + return d; +} + +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 +//============================================================================== + +class SurfaceXPlane : public Surface { + double x0; +public: + SurfaceXPlane(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + const bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; +}; + +SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node) { + const char *coeffs = surf_node.attribute("coeffs").value(); + int stat = sscanf(coeffs, "%lf", &x0); + if (stat != 1) { + std::cout << "Something went wrong reading surface coeffs!" << std::endl; + } +} + +inline double SurfaceXPlane::evaluate(const double xyz[3]) const { + return axis_aligned_plane_evaluate<0>(xyz, x0); +} + +inline double SurfaceXPlane::distance(const double xyz[3], const double uvw[3], + bool coincident) const { + return axis_aligned_plane_distance<0>(xyz, uvw, coincident, x0); +} + +inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const { + axis_aligned_plane_normal<0, 1, 2>(xyz, uvw); +} + +//============================================================================== +// SurfaceYPlane +//============================================================================== + +class SurfaceYPlane : public Surface { + double y0; +public: + 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; +}; + +SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) { + const char *coeffs = surf_node.attribute("coeffs").value(); + int stat = sscanf(coeffs, "%lf", &y0); + if (stat != 1) { + std::cout << "Something went wrong reading surface coeffs!" << std::endl; + } +} + +inline double SurfaceYPlane::evaluate(const double xyz[3]) const { + return axis_aligned_plane_evaluate<1>(xyz, y0); +} + +inline double SurfaceYPlane::distance(const double xyz[3], const double uvw[3], + bool coincident) const { + return axis_aligned_plane_distance<1>(xyz, uvw, coincident, y0); +} + +inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const { + axis_aligned_plane_normal<1, 0, 2>(xyz, uvw); +} + +//============================================================================== +// SurfaceZPlane //============================================================================== class SurfaceZPlane : public Surface { double z0; public: SurfaceZPlane(pugi::xml_node surf_node); - double evaluate(double xyz[3]); - double distance(double xyz[3], double uvw[3], bool coincident); - void normal(double xyz[3], double uvw[3]); + 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; }; SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) { @@ -95,73 +191,49 @@ SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) { } } -double -SurfaceZPlane::evaluate(double xyz[3]) { - double f = xyz[2] - z0; - return f; +inline double SurfaceZPlane::evaluate(const double xyz[3]) const { + return axis_aligned_plane_evaluate<2>(xyz, z0); } -double -SurfaceZPlane::distance(double xyz[3], double uvw[3], bool coincident) { - double f = z0 - xyz[2]; - double d; - if (coincident or fabs(f) < FP_COINCIDENT or uvw[2] == 0.0) - { - d = INFTY; - } - else - { - d = f / uvw[2]; - if (d < 0.0) d = INFTY; - } - return d; +inline double SurfaceZPlane::distance(const double xyz[3], const double uvw[3], + bool coincident) const { + return axis_aligned_plane_distance<2>(xyz, uvw, coincident, z0); } -void -SurfaceZPlane::normal(double xyz[3], double uvw[3]) { - uvw[0] = 0.0; - uvw[1] = 0.0; - uvw[2] = 1.0; +inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const { + axis_aligned_plane_normal<2, 0, 1>(xyz, uvw); } +//============================================================================== +// SurfacePlane //============================================================================== -class SurfaceZCylinder : public Surface { - double x0, y0, r; -public: - SurfaceZCylinder(pugi::xml_node surf_node); - double evaluate(double xyz[3]); - double distance(double xyz[3], double uvw[3], bool coincident); - void normal(double xyz[3], double uvw[3]); -}; +// TODO -SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node) { - const char *coeffs = surf_node.attribute("coeffs").value(); - int stat = sscanf(coeffs, "%lf %lf %lf", &x0, &y0, &r); - if (stat != 3) { - std::cout << "Something went wrong reading surface coeffs!" << std::endl; - } +//============================================================================== +// Generic functions for X-, Y-, and Z-, cylinders +//============================================================================== + +template double +axis_aligned_cylinder_evaluate(const double xyz[3], 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; } -double -SurfaceZCylinder::evaluate(double xyz[3]) { - const double x = xyz[0] - x0; - const double y = xyz[1] - y0; - double f = x*x + y*y - r*r; - return f; -} - -double -SurfaceZCylinder::distance(double xyz[3], double uvw[3], bool coincident) { - double a = 1.0 - uvw[2]*uvw[2]; // u^2 + v^2 +template double +axis_aligned_cylinder_distance(const double xyz[3], const double uvw[3], + double offset1, double offset2, double radius, bool coincident) { + const double a = 1.0 - uvw[i1]*uvw[i1]; // u^2 + v^2 if (a == 0.0) return INFTY; - double x = xyz[0] - x0; - double y = xyz[1] - y0; - double k = x*uvw[0] + y*uvw[1]; - double c = x*x + y*y - r*r; - double quad = k*k - a*c; + 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 quad = k*k - a*c; if (quad < 0.0) { // No intersection with cylinder. @@ -190,11 +262,123 @@ SurfaceZCylinder::distance(double xyz[3], double uvw[3], bool coincident) { } } -void -SurfaceZCylinder::normal(double xyz[3], double uvw[3]) { - uvw[0] = 2.0 * (xyz[0] - x0); - uvw[1] = 2.0 * (xyz[1] - y0); - uvw[2] = 0.0; +template void +axis_aligned_cylinder_normal(const double xyz[3], double uvw[3], double offset1, + double offset2) { + uvw[i2] = 2.0 * (xyz[i2] - offset1); + uvw[i3] = 2.0 * (xyz[i3] - offset2); + uvw[i1] = 0.0; +} + +//============================================================================== +// SurfaceXCylinder +//============================================================================== + +class SurfaceXCylinder : public Surface { + double y0, z0, r; +public: + 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; +}; + +SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node) { + const char *coeffs = surf_node.attribute("coeffs").value(); + int stat = sscanf(coeffs, "%lf %lf %lf", &y0, &z0, &r); + if (stat != 3) { + std::cout << "Something went wrong reading surface coeffs!" << std::endl; + } +} + +inline double +SurfaceXCylinder::evaluate(const double xyz[3]) const { + return axis_aligned_cylinder_evaluate<1, 2>(xyz, y0, z0, r); +} + +inline double SurfaceXCylinder::distance(const double xyz[3], + const double uvw[3], bool coincident) const { + return axis_aligned_cylinder_distance<0, 1, 2>(xyz, uvw, coincident, y0, z0, + r); +} + +inline void SurfaceXCylinder::normal(const double xyz[3], double uvw[3]) const { + axis_aligned_cylinder_normal<0, 1, 2>(xyz, uvw, y0, z0); +} + +//============================================================================== +// SurfaceYCylinder +//============================================================================== + +class SurfaceYCylinder : public Surface { + double x0, z0, r; +public: + 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; +}; + +SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node) { + const char *coeffs = surf_node.attribute("coeffs").value(); + int stat = sscanf(coeffs, "%lf %lf %lf", &x0, &z0, &r); + if (stat != 3) { + std::cout << "Something went wrong reading surface coeffs!" << std::endl; + } +} + +inline double +SurfaceYCylinder::evaluate(const double xyz[3]) const { + return axis_aligned_cylinder_evaluate<0, 2>(xyz, x0, z0, r); +} + +inline double SurfaceYCylinder::distance(const double xyz[3], + const double uvw[3], bool coincident) const { + return axis_aligned_cylinder_distance<1, 0, 2>(xyz, uvw, coincident, x0, z0, + r); +} + +inline void SurfaceYCylinder::normal(const double xyz[3], double uvw[3]) const { + axis_aligned_cylinder_normal<1, 0, 2>(xyz, uvw, x0, z0); +} + +//============================================================================== +// SurfaceZCylinder +//============================================================================== + +class SurfaceZCylinder : public Surface { + double x0, y0, r; +public: + 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; +}; + +SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node) { + const char *coeffs = surf_node.attribute("coeffs").value(); + int stat = sscanf(coeffs, "%lf %lf %lf", &x0, &y0, &r); + if (stat != 3) { + std::cout << "Something went wrong reading surface coeffs!" << std::endl; + } +} + +inline double +SurfaceZCylinder::evaluate(const double xyz[3]) const { + return axis_aligned_cylinder_evaluate<0, 1>(xyz, x0, y0, r); +} + +inline double SurfaceZCylinder::distance(const double xyz[3], + const double uvw[3], bool coincident) const { + return axis_aligned_cylinder_distance<2, 0, 1>(xyz, uvw, coincident, x0, y0, + r); +} + +inline void SurfaceZCylinder::normal(const double xyz[3], double uvw[3]) const { + axis_aligned_cylinder_normal<2, 0, 1>(xyz, uvw, x0, y0); } //============================================================================== @@ -219,16 +403,26 @@ read_surfaces(pugi::xml_node *node) { surf_node = surf_node.next_sibling("surface"), i_surf++) { if (surf_node.attribute("type")) { const pugi::char_t *surf_type = surf_node.attribute("type").value(); - if (strcmp(surf_type, "z-cylinder") == 0) - { - surfaces_c[i_surf] = new SurfaceZCylinder(surf_node); - } - else if (strcmp(surf_type, "z-plane") == 0) - { + + if (strcmp(surf_type, "x-plane") == 0) { + surfaces_c[i_surf] = new SurfaceXPlane(surf_node); + + } else if (strcmp(surf_type, "y-plane") == 0) { + surfaces_c[i_surf] = new SurfaceYPlane(surf_node); + + } else if (strcmp(surf_type, "z-plane") == 0) { surfaces_c[i_surf] = new SurfaceZPlane(surf_node); - } - else - { + + } else if (strcmp(surf_type, "x-cylinder") == 0) { + surfaces_c[i_surf] = new SurfaceXCylinder(surf_node); + + } else if (strcmp(surf_type, "y-cylinder") == 0) { + surfaces_c[i_surf] = new SurfaceYCylinder(surf_node); + + } else if (strcmp(surf_type, "z-cylinder") == 0) { + surfaces_c[i_surf] = new SurfaceZCylinder(surf_node); + + } else { std::cout << "Call error or handle uppercase here!" << std::endl; std::cout << surf_type << std::endl; }