geometry/shapes/cylinder3D.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 // Project include(s)
0012 #include "detray/definitions/algebra.hpp"
0013 #include "detray/definitions/containers.hpp"
0014 #include "detray/definitions/detail/qualifiers.hpp"
0015 #include "detray/definitions/indexing.hpp"
0016 #include "detray/definitions/math.hpp"
0017 #include "detray/definitions/units.hpp"
0018 #include "detray/geometry/coordinates/cylindrical3D.hpp"
0019 #include "detray/geometry/detail/shape_utils.hpp"
0020
0021 // System include(s)
0022 #include <limits>
0023 #include <ostream>
0024 #include <string_view>
0025
0026 namespace detray {
0027
0028 /// @brief Geometrical shape of a full 3D cylinder.
0029 ///
0030 /// It is defined by r and the two half lengths rel to the coordinate center.
0031 class cylinder3D {
0032  public:
0033   /// The name for this shape
0034   static constexpr std::string_view name = "cylinder3D";
0035
0036   enum boundaries : unsigned int {
0037     e_min_r = 0u,
0038     e_min_phi = 1u,
0039     e_min_z = 2u,
0040     e_max_r = 3u,
0041     e_max_phi = 4u,
0042     e_max_z = 5u,
0043     e_size = 6u,
0044   };
0045
0046   /// Container definition for the shape boundary values
0047   template <concepts::scalar scalar_t>
0048   using bounds_type = darray<scalar_t, boundaries::e_size>;
0049
0050   /// Local coordinate frame for boundary checks
0051   template <concepts::algebra algebra_t>
0052   using local_frame_type = cylindrical3D<algebra_t>;
0053
0054   /// Result type of a boundary check
0055   template <typename bool_t>
0056   using result_type = bool_t;
0057
0058   /// Dimension of the local coordinate system
0059   static constexpr std::size_t dim{3u};
0060
0061   /// @brief Find the minimum distance to any boundary.
0062   ///
0063   /// @note the point is expected to be given in local coordinates by the
0064   /// caller.
0065   ///
0066   /// @param bounds the boundary values for this shape
0067   /// @param loc_p the point to be checked in the local coordinate system
0068   ///
0069   /// @return the minimum distance.
0070   template <concepts::scalar scalar_t, concepts::point point_t>
0071   DETRAY_HOST_DEVICE inline scalar_t min_dist_to_boundary(
0072       const bounds_type<scalar_t> &bounds, const point_t &loc_p) const {
0073     const scalar_t min_r_dist =
0074         math::min(math::fabs(loc_p[0] - bounds[e_min_r]),
0075                   math::fabs(bounds[e_max_r] - loc_p[0]));
0076     const scalar_t min_phi_dist =
0077         math::min(math::fabs(loc_p[1] - bounds[e_min_phi]),
0078                   math::fabs(bounds[e_max_phi] - loc_p[1]));
0079     const scalar_t min_z_dist =
0080         math::min(math::fabs(loc_p[2] - bounds[e_min_z]),
0081                   math::fabs(bounds[e_max_z] - loc_p[2]));
0082
0083     // Use the chord for the phi distance
0084     return math::min(
0085         math::min(min_r_dist, 2.f * loc_p[0] * math::sin(0.5f * min_phi_dist)),
0086         min_z_dist);
0087   }
0088
0089   /// @brief Check boundary values for a local point.
0090   /// @{
0091   /// @param bounds the boundary values for this shape
0092   /// @param trf the surface transform
0093   /// @param glob_p the point to be checked in the global coordinate system
0094   /// @param tol dynamic tolerance determined by caller
0095   ///
0096   /// @return true if the local point lies within the given boundaries.
0097   template <concepts::algebra algebra_t>
0098   DETRAY_HOST_DEVICE constexpr result_type<dbool<algebra_t>> check_boundaries(
0099       const bounds_type<dscalar<algebra_t>> &bounds,
0100       const dtransform3D<algebra_t> &trf, const dpoint3D<algebra_t> &glob_p,
0101       const dscalar<algebra_t> tol =
0102           std::numeric_limits<dscalar<algebra_t>>::epsilon(),
0103       const dscalar<algebra_t> /*edge_tol*/ = 0.f) const {
0104     // Get the full local position
0105     const auto loc_p =
0106         local_frame_type<algebra_t>::global_to_local(trf, glob_p, {});
0107
0108     return check_boundaries(bounds, loc_p, tol);
0109   }
0110
0111   /// @note the point is expected to be given in local coordinates by the
0112   /// caller. For the conversion from global cartesian coordinates, the
0113   /// nested @c shape struct can be used. The point is assumed to be in
0114   /// the cylinder 3D frame.
0115   ///
0116   /// @param bounds the boundary values for this shape
0117   /// @param loc_p the point to be checked in the local coordinate system
0118   /// @param tol dynamic tolerance determined by caller
0119   ///
0120   /// @return true if the local point lies within the given boundaries.
0121   template <concepts::scalar scalar_t, concepts::point point_t>
0122   DETRAY_HOST_DEVICE constexpr auto check_boundaries(
0123       const bounds_type<scalar_t> &bounds, const point_t &loc_p,
0124       const scalar_t tol = std::numeric_limits<scalar_t>::epsilon(),
0125       const scalar_t /*edge_tol*/ = 0.f) const {
0126     const scalar_t phi_tol = detail::phi_tolerance(tol, loc_p[0]);
0127
0128     return ((bounds[e_min_r] - tol) <= loc_p[0] &&
0129             (bounds[e_min_phi] - phi_tol) <= loc_p[1] &&
0130             (bounds[e_min_z] - tol) <= loc_p[2] &&
0131             loc_p[0] <= (bounds[e_max_r] + tol) &&
0132             loc_p[1] <= (bounds[e_max_phi] + phi_tol) &&
0133             loc_p[2] <= (bounds[e_max_z] + tol));
0134   }
0135   /// @}
0136
0137   /// @brief Measure of the shape: Volume
0138   ///
0139   /// @param bounds the boundary values for this shape
0140   ///
0141   /// @returns the cylinder volume as parto of global space.
0142   template <concepts::scalar scalar_t>
0143   DETRAY_HOST_DEVICE constexpr scalar_t measure(
0144       const bounds_type<scalar_t> &bounds) const {
0145     return volume(bounds);
0146   }
0147
0148   /// @brief The volume of a the shape
0149   ///
0150   /// @param bounds the boundary values for this shape
0151   ///
0152   /// @returns the cylinder volume.
0153   template <concepts::scalar scalar_t>
0154   DETRAY_HOST_DEVICE constexpr scalar_t volume(
0155       const bounds_type<scalar_t> &bounds) const {
0156     return constant<scalar_t>::pi * (bounds[e_max_z] - bounds[e_min_z]) *
0157            (bounds[e_max_r] * bounds[e_max_r] -
0158             bounds[e_min_r] * bounds[e_min_r]);
0159   }
0160
0161   /// @brief Merge two 3D cylinder shapes
0162   ///
0163   /// @param bounds the boundary values for this shape
0164   /// @param o_bounds the boundary values for the other shape
0165   ///
0166   /// @returns merged bound values
0167   template <concepts::scalar scalar_t>
0168   DETRAY_HOST_DEVICE constexpr bounds_type<scalar_t> merge(
0169       const bounds_type<scalar_t> &bounds,
0170       const bounds_type<scalar_t> &o_bounds) const {
0171     bounds_type<scalar_t> new_bounds{};
0172
0173     new_bounds[e_min_r] = math::min(bounds[e_min_r], o_bounds[e_min_r]);
0174     new_bounds[e_min_phi] = math::min(bounds[e_min_phi], o_bounds[e_min_phi]);
0175     new_bounds[e_min_z] = math::min(bounds[e_min_z], o_bounds[e_min_z]);
0176     new_bounds[e_max_r] = math::max(bounds[e_max_r], o_bounds[e_max_r]);
0177     new_bounds[e_max_phi] = math::max(bounds[e_max_phi], o_bounds[e_max_phi]);
0178     new_bounds[e_max_z] = math::max(bounds[e_max_z], o_bounds[e_max_z]);
0179
0180     return new_bounds;
0181   }
0182
0183   /// @brief Lower and upper point for minimal axis aligned bounding box.
0184   ///
0185   /// Computes the min and max vertices in a local cartesian frame.
0186   ///
0187   /// @param bounds the boundary values for this shape
0188   /// @param env dynamic envelope around the shape
0189   ///
0190   /// @returns and array of coordinates that contains the lower point (first
0191   /// three values) and the upper point (latter three values).
0192   // @todo: Look at phi - range for a better fit
0193   template <concepts::algebra algebra_t>
0194   DETRAY_HOST_DEVICE inline darray<dscalar<algebra_t>, 6> local_min_bounds(
0195       const bounds_type<dscalar<algebra_t>> &bounds,
0196       const dscalar<algebra_t> env =
0197           std::numeric_limits<dscalar<algebra_t>>::epsilon()) const {
0198     assert(env > 0.f);
0199     const dscalar<algebra_t> r_bound{bounds[e_max_r] + env};
0200     return {-r_bound, -r_bound, bounds[e_min_z] - env,
0201             r_bound,  r_bound,  bounds[e_max_z] + env};
0202   }
0203
0204   /// @returns the shapes centroid in local cartesian coordinates
0205   template <concepts::algebra algebra_t>
0206   DETRAY_HOST_DEVICE dpoint3D<algebra_t> centroid(
0207       const bounds_type<dscalar<algebra_t>> &bounds) const {
0208     return 0.5f * dpoint3D<algebra_t>{static_cast<dscalar<algebra_t>>(0.f),
0209                                       (bounds[e_min_phi] + bounds[e_max_phi]),
0210                                       (bounds[e_min_z] + bounds[e_max_z])};
0211   }
0212
0213   /// Generate vertices in local cartesian frame
0214   ///
0215   /// @param bounds the boundary values for the cylinder
0216   /// @param n_seg is the number of line segments
0217   ///
0218   /// @return a generated list of vertices
0219   template <concepts::algebra algebra_t>
0220   DETRAY_HOST dvector<dpoint3D<algebra_t>> vertices(
0221       const bounds_type<dscalar<algebra_t>> & /*unused*/,
0222       dindex /*unused*/) const {
0223     throw std::runtime_error(
0224         "Vertex generation for 3D cylinders is not implemented");
0225     return {};
0226   }
0227
0228   /// @brief Check consistency of boundary values.
0229   ///
0230   /// @param bounds the boundary values for this shape
0231   /// @param os output stream for error messages
0232   ///
0233   /// @return true if the bounds are consistent.
0234   template <concepts::scalar scalar_t>
0235   DETRAY_HOST constexpr bool check_consistency(
0236       const bounds_type<scalar_t> &bounds, std::ostream &os) const {
0237     constexpr auto tol{10.f * std::numeric_limits<scalar_t>::epsilon()};
0238
0239     if (bounds[e_min_r] < tol) {
0240       os << "DETRAY ERROR (HOST): Radii must be in the range (0, "
0241             "numeric_max)"
0242          << std::endl;
0243       return false;
0244     }
0245     if (bounds[e_min_r] >= bounds[e_max_r] ||
0246         math::fabs(bounds[e_min_r] - bounds[e_max_r]) < tol) {
0247       os << "DETRAY ERROR (HOST): Min Radius must be smaller than max "
0248             "Radius.";
0249       return false;
0250     }
0251     if (bounds[e_min_z] >= bounds[e_max_z] ||
0252         math::fabs(bounds[e_min_z] - bounds[e_max_z]) < tol) {
0253       os << "DETRAY ERROR (HOST): Min z must be smaller than max z.";
0254       return false;
0255     }
0256
0257     return true;
0258   }
0259 };
0260
0261 }  // namespace detray