orange/surf/Involute.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file orange/surf/Involute.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cassert>
0010 #include <cmath>
0011
0012 #include "corecel/Constants.hh"
0013 #include "corecel/Types.hh"
0014 #include "corecel/cont/Array.hh"
0015 #include "corecel/cont/Span.hh"
0016 #include "corecel/math/Algorithms.hh"
0017 #include "corecel/math/ArrayOperators.hh"
0018 #include "corecel/math/ArrayUtils.hh"
0019 #include "orange/OrangeTypes.hh"
0020
0021 #include "detail/InvolutePoint.hh"
0022 #include "detail/InvoluteSolver.hh"
0023
0024 namespace celeritas
0025 {
0026 //---------------------------------------------------------------------------//
0027 /*!
0028  * Z-aligned circular involute.
0029  *
0030  * An involute is a curve created by unwinding a chord from a shape.
0031  * The involute implemented here is constructed from a circle and can be made
0032  * to be clockwise (negative) or anti-clockwise (positive).
0033  * This involute is the same type of that found in HFIR, and is necessary for
0034  * building accurate models.
0035  *
0036  * While the parameters define the curve itself have no restriction, except for
0037  * the radius of involute being positive, the involute needs to be bounded for
0038  * there to be finite solutions to the roots of the intersection points.
0039  * Additionally this bound cannot exceed an interval of size of 2 to be able to
0040  * determine if whether a particle is in, out, or on the involute surface.
0041  * Lastly the involute geometry is fixed to be z-aligned.
0042  *
0043  * \f[
0044  *   x = r_b (\cos(t+a) + t \sin(t+a)) + x_0
0045  * \f]
0046  * \f[
0047  *   y = r_b (\sin(t+a) - t \cos(t+a)) + y_0
0048  * \f]
0049  *
0050  * where \em t is the normal angle of the tangent to the circle of involute
0051  * with radius \em r_b from a starting angle of \em a (\f$r/h\f$ for a finite
0052  * cone.
0053  */
0054 class Involute
0055 {
0056   public:
0057     //@{
0058     //! \name Type aliases
0059     using Intersections = Array<real_type, 3>;
0060     using StorageSpan = Span<real_type const, 6>;
0061     using Real2 = Array<real_type, 2>;
0062     //@}
0063
0064   public:
0065     //// CLASS ATTRIBUTES ////
0066
0067     // Surface type identifier
0068     static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type()
0069
0070     {
0071         return SurfaceType::inv;
0072     }
0073
0074     //! Safety cannot be trivially calculated
0075     static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return false; }
0076
0077   public:
0078     //// CONSTRUCTORS ////
0079
0080     explicit Involute(Real2 const& origin,
0081                       real_type radius,
0082                       real_type displacement,
0083                       Chirality sign,
0084                       real_type tmin,
0085                       real_type tmax);
0086
0087     // Construct from raw data
0088     template<class R>
0089     explicit inline CELER_FUNCTION Involute(Span<R, StorageSpan::extent>);
0090
0091     //// ACCESSORS ////
0092
0093     //! X-Y center of the circular base of the involute
0094     CELER_FUNCTION Real2 const& origin() const { return origin_; }
0095
0096     //! Involute circle's radius
0097     CELER_FUNCTION real_type r_b() const { return std::fabs(r_b_); }
0098
0099     //! Displacement angle
0100     CELER_FUNCTION real_type displacement_angle() const { return a_; }
0101
0102     // Orientation of the involute curve
0103     inline CELER_FUNCTION Chirality sign() const;
0104
0105     //! Get bounds of the involute
0106     CELER_FUNCTION real_type tmin() const { return tmin_; }
0107     CELER_FUNCTION real_type tmax() const { return tmax_; }
0108
0109     //! Get a view to the data for type-deleted storage
0110     CELER_FUNCTION StorageSpan data() const { return {&origin_[0], 6}; }
0111
0112     //// CALCULATION ////
0113
0114     // Determine the sense of the position relative to this surface
0115     inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const;
0116
0117     // Calculate all possible straight-line intersections with this surface
0118     inline CELER_FUNCTION Intersections calc_intersections(
0119         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0120
0121     // Calculate outward normal at a position
0122     inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const;
0123
0124   private:
0125     // Location of the center of circle of involute
0126     Real2 origin_;
0127
0128     // Involute parameters
0129     real_type r_b_;  // Radius, negative if "clockwise" (flipped)
0130     real_type a_;
0131
0132     // Bounds
0133     real_type tmin_;
0134     real_type tmax_;
0135 };
0136
0137 //---------------------------------------------------------------------------//
0138 // INLINE DEFINITIONS
0139 //---------------------------------------------------------------------------//
0140 /*!
0141  * Construct from raw data.
0142  */
0143 template<class R>
0144 CELER_FUNCTION Involute::Involute(Span<R, StorageSpan::extent> data)
0145     : origin_{data[0], data[1]}
0146     , r_b_{data[2]}
0147     , a_{data[3]}
0148     , tmin_{data[4]}
0149     , tmax_{data[5]}
0150 {
0151 }
0152
0153 //---------------------------------------------------------------------------//
0154 /*!
0155  * Orientation of the involute curve.
0156  */
0157 CELER_FUNCTION auto Involute::sign() const -> Chirality
0158 {
0159     return r_b_ > 0 ? Chirality::left : Chirality::right;
0160 }
0161
0162 //---------------------------------------------------------------------------//
0163 /*!
0164  * Determine the sense of the position relative to this surface.
0165  *
0166  * This is performed by first checking if the particle is outside the bounding
0167  * circles given by tmin and tmax.
0168  *
0169  * Then position is compared to a point on the involute at the same radius from
0170  * the origin.
0171  *
0172  * Next the point is determined to be inside if an involute that passes it,
0173  * has an \em t that is within the bounds and an \em a greater than the
0174  * involute being tested for. To do this, the tangent to the involute must be
0175  * determined by calculating the intercepts of the circle of involute and
0176  * circle centered at the midpoint of the origin and the point being tested for
0177  * that extends to the origin. After which the angle of tangent point can be
0178  * determined and used to calculate \em a via: \f[
0179  * a = angle - t_point
0180  * \f]
0181  *
0182  * To calculate the intercept of two circles it is assumed that the two circles
0183  * are \em x_prime axis. The coordinate along \em x_prime originating from the
0184  * center of one of the circles is given by: \f[ x_prime = (d^2 - r^2 +
0185  * R^2)/(2d) \f] Where \em d is the distance between the center of the two
0186  * circles, \em r is the radius of the displaced circle, and \em R is the
0187  * radius of the circle at the orgin. Then the point \em y_prime can be
0188  * obatined from: \f[ y_prime = \pm sqrt(R^2 - x_prime^2) \f] Then to convert (
0189  * \em x_prime , \em y_prime ) to ( \em x, \em y ) the following rotation is
0190  * applied:\f[ x = x_prime*cos(theta) - y_prime*sin(theta) y =
0191  * y_prime*cos(theta) + x_prime*sin(theta) \f] Where \em theta is the angle
0192  * between the \em x_prime axis and the \em x axis.
0193  */
0194 CELER_FUNCTION SignedSense Involute::calc_sense(Real3 const& pos) const
0195 {
0196     using constants::pi;
0197
0198     Real2 xy;
0199     xy[0] = pos[0] - origin_[0];
0200     xy[1] = pos[1] - origin_[1];
0201
0202     if (this->sign() == Chirality::right)
0203     {
0204         xy[0] = negate(xy[0]);
0205     }
0206
0207     // Calculate distance to origin and obtain t value for distance.
0208     real_type const t_point_sq = dot_product(xy, xy) / ipow<2>(r_b_) - 1;
0209
0210     // Check if point is in defined bounds.
0211     if (t_point_sq < ipow<2>(tmin_))
0212     {
0213         return SignedSense::outside;
0214     }
0215     if (t_point_sq > ipow<2>(tmax_))
0216     {
0217         return SignedSense::outside;
0218     }
0219
0220     // Check if Point is on involute.
0221     detail::InvolutePoint calc_point{this->r_b(), a_};
0222     Real2 point = calc_point(clamp_to_nonneg(t_point_sq));
0223
0224     if (xy == point)
0225     {
0226         return SignedSense::on;
0227     }
0228
0229     // Check if point is in interval
0230
0231     // Calculate tangent point
0232     // TODO: check that compiler avoids recomputing trig functions
0233     real_type x_prime = ipow<2>(r_b_) / norm(xy);
0234     real_type y_prime = std::sqrt(ipow<2>(r_b_) - ipow<2>(x_prime));
0235
0236     // TODO: check that compiler avoids recomputing trig functions
0237     point[0] = (x_prime * xy[0] - y_prime * xy[1]) / norm(xy);
0238     point[1] = (y_prime * xy[0] + x_prime * xy[1]) / norm(xy);
0239
0240     // Calculate angle of tangent
0241     real_type theta = std::acos(point[0] / norm(point));
0242     if (point[1] < 0)
0243     {
0244         theta = 2 * pi - theta;
0245     }
0246     // Count number of positive rotations around involute
0247     theta += max<real_type>(real_type{0},
0248                             std::floor((tmax_ + a_ - theta) / (2 * pi)))
0249              * 2 * pi;
0250
0251     // Calculate the displacement angle of the point
0252     real_type a1 = theta - std::sqrt(clamp_to_nonneg(t_point_sq));
0253
0254     // Check if point is inside bounds
0255     if (theta < tmax_ + a_ && a1 > a_)
0256     {
0257         return SignedSense::inside;
0258     }
0259
0260     return SignedSense::outside;
0261 }
0262
0263 //---------------------------------------------------------------------------//
0264 /*!
0265  * Calculate all possible straight-line intersections with this surface.
0266  */
0267 CELER_FUNCTION auto
0268 Involute::calc_intersections(Real3 const& pos,
0269                              Real3 const& dir,
0270                              SurfaceState on_surface) const -> Intersections
0271 {
0272     // Expand translated positions into 'xyz' coordinate system
0273     Real3 rel_pos{pos};
0274     rel_pos[0] -= origin_[0];
0275     rel_pos[1] -= origin_[1];
0276
0277     detail::InvoluteSolver solve(this->r_b(), a_, this->sign(), tmin_, tmax_);
0278
0279     return solve(rel_pos, dir, on_surface);
0280 }
0281
0282 //---------------------------------------------------------------------------//
0283 /*!
0284  * Calculate outward normal at a position on the surface.
0285  *
0286  * This is done by taking the derivative of the involute equation : \f[
0287  * n = {\sin(t+a) - \cos(t+a), 0}
0288  * \f]
0289  */
0290 CELER_FORCEINLINE_FUNCTION Real3 Involute::calc_normal(Real3 const& pos) const
0291 {
0292     Real2 xy;
0293     xy[0] = pos[0] - origin_[0];
0294     xy[1] = pos[1] - origin_[1];
0295
0296     // Calculate normal
0297     real_type const angle
0298         = std::sqrt(clamp_to_nonneg(dot_product(xy, xy) / ipow<2>(r_b_) - 1))
0299           + a_;
0300     Real3 normal_ = {std::sin(angle), -std::cos(angle), 0};
0301
0302     if (this->sign() == Chirality::right)
0303     {
0304         normal_[0] = negate(normal_[0]);
0305     }
0306
0307     return normal_;
0308 }
0309
0310 //---------------------------------------------------------------------------//
0311 }  // namespace celeritas