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