Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <limits>
       9             : #include <optional>
      10             : 
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : #include "Utilities/Gsl.hpp"
      13             : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
      14             : 
      15             : /// Holds helper functions for use with
      16             : /// domain::CoordinateMaps::FocallyLiftedMap.
      17           1 : namespace domain::CoordinateMaps::FocallyLiftedMapHelpers {
      18             : 
      19             : /*!
      20             :  * \brief Finds how long to extend a line segment to have it intersect
      21             :  * a point on a 2-sphere.
      22             :  *
      23             :  * \details Consider a 2-sphere with center \f$C^i\f$ and radius \f$R\f$, and
      24             :  * and let \f$P^i\f$ and \f$x_0^i\f$ be two arbitrary 3D points.
      25             :  *
      26             :  * Consider the line passing through \f$P^i\f$ and \f$x_0^i\f$.
      27             :  * If this line intersects the sphere at a point \f$x_1^i\f$, then we can write
      28             :  *
      29             :  * \f[
      30             :  *  x_1^i = P^i + (x_0^i-P^i) \lambda,
      31             :  * \f]
      32             :  *
      33             :  * where \f$\lambda\f$ is a scale factor.
      34             :  *
      35             :  * `scale_factor` computes and returns \f$\lambda\f$.
      36             :  *
      37             :  * ### Even more detail:
      38             :  *
      39             :  * To solve for \f$\lambda\f$, we note that \f$x_1^i\f$ is on the surface of
      40             :  * the sphere, so
      41             :  *
      42             :  * \f[
      43             :  *  |x_1^i-C^i|^2 = R^2,
      44             :  * \f]
      45             :  *
      46             :  * (where \f$|A^i|^2\f$ means \f$\delta_{ij} A^i A^j\f$),
      47             :  *
      48             :  *  or equivalently
      49             :  *
      50             :  * \f[
      51             :  *  | P^i-C^i + (x_0^i-P^i)\lambda |^2 = R^2.
      52             :  * \f]
      53             :  *
      54             :  * This is a quadratic equation for \f$\lambda\f$
      55             :  * and it generally has more than one real root.
      56             :  * It takes the usual form \f$a\lambda^2+b\lambda+c=0\f$,
      57             :  * with
      58             :  *
      59             :  * \f{align*}
      60             :  *  a &= |x_0^i-P^i|^2,\\
      61             :  *  b &= 2(x_0^i-P^i)(P^j-C^j)\delta_{ij},\\
      62             :  *  c &= |P^i-C^i|^2 - R^2,
      63             :  * \f}
      64             :  *
      65             :  * So how do we choose between multiple roots?  Some of the maps that
      66             :  * use `scale_factor` assume that *for all points*, \f$x_0^i\f$ is
      67             :  * between \f$P^i\f$ and \f$x_1^i\f$.  Those maps should set the parameter
      68             :  * `src_is_between_proj_and_target` to true. Other maps assume that
      69             :  * *for all points*, \f$x^i\f$ is always between \f$x_0^i\f$
      70             :  * and \f$P^i\f$. Those maps should set the parameter
      71             :  * `src_is_between_proj_and_target` to false.
      72             :  *
      73             :  * \warning If we ever add maps where
      74             :  * `src_is_between_proj_and_target` can change from point to point,
      75             :  * the logic of `scale_factor` needs to be changed.
      76             :  *
      77             :  * In the arguments to the function below, `src_point` is  \f$x_0^i\f$,
      78             :  * `proj_center` is \f$P^i\f$, `sphere_center` is \f$C^i\f$, and
      79             :  * `radius` is \f$R\f$.
      80             :  *
      81             :  */
      82             : template <typename T>
      83           1 : void scale_factor(const gsl::not_null<tt::remove_cvref_wrap_t<T>*>& result,
      84             :                   const std::array<T, 3>& src_point,
      85             :                   const std::array<double, 3>& proj_center,
      86             :                   const std::array<double, 3>& sphere_center, double radius,
      87             :                   bool src_is_between_proj_and_target);
      88             : 
      89             : /*!
      90             :  *  Solves a problem of the same form as `scale_factor`, but is used
      91             :  *  only by the inverse function to compute \f$\tilde{\lambda}\f$ and
      92             :  *  \f$\bar{\lambda}\f$. `try_scale_factor` is used in two contexts:
      93             :  *
      94             :  *  `try_scale_factor` is used to determine \f$\bar{\lambda}\f$
      95             :  *  given \f$x^i\f$. \f$\bar{\lambda}\f$ is defined by
      96             :  *  \f{align*} x_1^i = P^i + (x^i - P^i) \bar{\lambda}.\f}
      97             :  *
      98             :  *  `try_scale_factor` is used by the `lambda_tilde` functions of some
      99             :  *  of the `InnerMap` classes (namely those `InnerMap` classes where
     100             :  *  \f$x_0^i\f$ is a spherical surface) to solve for
     101             :  *  \f$\tilde{\lambda}\f$ given\f$x^i\f$.  \f$\tilde{\lambda}\f$
     102             :  *  is defined by
     103             :  *  \f{align*} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\f}
     104             :  *
     105             :  *  In both of these contexts, the input parameter `src_point` is
     106             :  *  \f$x^i\f$, a point that is supposed to be in the range of the
     107             :  *  `FocallyLiftedMap`. Because the inverse function can be and is
     108             :  *  called for an arbitrary \f$x^i\f$ that might not be in the range
     109             :  *  of the `FocallyLiftedMap`, `try_scale_factor` returns a
     110             :  *  std::optional, with a default-constructed std::optional if the roots it
     111             :  *  finds are not as expected (i.e. if the inverse map was called for
     112             :  *  a point not in the range of the map).
     113             :  *
     114             :  *  Because `try_scale_factor` can be called in different situations,
     115             :  *  it has additional boolean arguments `pick_larger_root` and
     116             :  *  `pick_root_greater_than_one` that allow the caller to choose which
     117             :  *  root to return.
     118             :  *
     119             :  *  Furthermore, to reduce roundoff errors near
     120             :  *  \f$\tilde{\lambda}=1\f$, the default behavior is to solve the
     121             :  *  quadratic equation for \f$\tilde{\lambda}-1\f$ (and then add
     122             :  *  \f$1\f$ to the solution). If instead one wants to solve the
     123             :  *  quadratic equation directly for \f$\tilde{\lambda}\f$ so as to
     124             :  *  obtain slightly different roundoff behavior, then one should
     125             :  *  specify the argument `solve_for_root_minus_one` to be `false`.
     126             :  *
     127             :  * `try_scale_factor` is not templated
     128             :  *  on type because it is used only by the inverse function, which
     129             :  *  works only on doubles.
     130             :  *
     131             :  */
     132           1 : std::optional<double> try_scale_factor(
     133             :     const std::array<double, 3>& src_point,
     134             :     const std::array<double, 3>& proj_center,
     135             :     const std::array<double, 3>& sphere_center, double radius,
     136             :     bool pick_larger_root, bool pick_root_greater_than_one,
     137             :     bool solve_for_root_minus_one = true);
     138             : 
     139             : /*!
     140             :  * Computes \f$\partial \lambda/\partial x_0^i\f$, where \f$\lambda\f$
     141             :  * is the quantity returned by `scale_factor` and `x_0` is `src_point` in
     142             :  * the `scale_factor` function.
     143             :  *
     144             :  * The formula (see `FocallyLiftedMap`) is
     145             :  * \f{align*}
     146             :  * \frac{\partial\lambda}{\partial x_0^j} &=
     147             :  * \lambda^2 \frac{C_j - x_1^j}{|x_1^i - P^i|^2
     148             :  * + (x_1^i - P^i)(P_i - C_i)}.
     149             :  * \f}
     150             :  *
     151             :  * Note that it takes `intersection_point` and not `src_point` as a
     152             :  * parameter.
     153             :  *
     154             :  * In the arguments to the function below, `intersection_point` is \f$x_1\f$,
     155             :  * `proj_center` is \f$P^i\f$, `sphere_center` is \f$C^i\f$, and
     156             :  * `radius` is \f$R\f$.
     157             :  */
     158             : template <typename T>
     159           1 : void d_scale_factor_d_src_point(
     160             :     const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>& result,
     161             :     const std::array<T, 3>& intersection_point,
     162             :     const std::array<double, 3>& proj_center,
     163             :     const std::array<double, 3>& sphere_center, const T& lambda);
     164             : 
     165             : }  // namespace domain::CoordinateMaps::FocallyLiftedMapHelpers