Line data    Source code

       1           0 : \cond NEVER
       2             : Distributed under the MIT License.
       3             : See LICENSE.txt for details.
       4             : \endcond
       5             : # Redistributing Gridpoints {#redistributing_gridpoints}
       6             : 
       7             : \tableofcontents
       8             : 
       9             : ## Introduction
      10             : The simplest way to construct a volume map from two parameterized surfaces is
      11             : by linearly interpolating between them:
      12             : 
      13             : \f[\vec{x}(\xi,\eta,\zeta) =
      14             : \frac{1-\zeta}{2}\vec{\sigma}_-(\xi,\eta)+
      15             : \frac{1+\zeta}{2}\vec{\sigma}_+(\xi,\eta)\f]
      16             : 
      17             : In the above example, each surface \f$\vec{\sigma}_+\f$ and
      18             : \f$\vec{\sigma}_-\f$ is parameterized using the logical coordinates \f$\xi\f$
      19             : and \f$\eta\f$, and a third coordinate \f$\zeta\in[-1,1]\f$ is used to
      20             : interpolate between them.
      21             : 
      22             : We then distribute gridpoints on this volume by specifying values of the
      23             : coordinates \f$\xi,\eta,\f$ and \f$\zeta\f$ at which the gridpoints are located.
      24             : In SpECTRE these values are the locations of the quadrature nodes. The
      25             : distribution of the gridpoints throughout the volume depends on the
      26             : parameterization used, and the simplest choice of parameterization does not
      27             : necessarily lead to the best gridpoint distribution. In this section we discuss
      28             : situations in which there exist better parameterizations than those obtained by
      29             : linear interpolation.
      30             : 
      31             : ## Generalized Logical Coordinates
      32             : 
      33             : In each of the following examples, we will obtain functions \f$\Xi(\xi),
      34             : \mathrm{H}(\eta),\f$ and \f$\mathrm{Z}(\zeta)\f$ that give better gridpoint
      35             : distributions than using the logical coordinates alone. Where possible, we will
      36             : write the reparameterized map such that the functional form of the map is
      37             : unchanged when replacing \f$\Xi\f$ with \f$\xi\f$, etc. We therefore refer to
      38             : \f$\Xi, \mathrm{H},\f$ and \f$\mathrm{Z}\f$ as the
      39             : *generalized logical coordinates*, as they can also refer to the logical
      40             : coordinates \f$\xi, \eta,\f$ and \f$\zeta\f$ themselves, when the transformation
      41             : is the identity.
      42             : 
      43             : ## Equiangular Maps
      44             : 
      45             : The mapping for a cubed sphere surface can be easily obtained by taking points
      46             : that lie on each face of a cube and normalizing them such that they lie on the
      47             : sphere:
      48             : 
      49             : \f[\vec{\sigma}_{+z}(\xi,\eta) =
      50             : \frac{1}{\sqrt{1 + \xi^2  + \eta^2}}
      51             : \begin{bmatrix}
      52             : \xi\\
      53             : \eta\\
      54             : 1\\
      55             : \end{bmatrix}\f]
      56             : 
      57             : In the above example the parameterization used for the upper \f$+z\f$ surface
      58             : of the cube is linear in \f$\xi\f$ and \f$\eta\f$. However, distances measured
      59             : on the surface of the sphere are not linear in \f$\xi\f$ and \f$\eta\f$. To see
      60             : this, one may compute \f$g_{\xi\xi} = |\frac{\partial\vec{x}}{\partial\xi}|^2\f$
      61             : to see how distances are measured in terms of \f$\xi\f$:
      62             : 
      63             : \f[g_{\xi,\xi}|_{\eta=0} = \frac{1}{(1+\xi^2)^2}\f]
      64             : 
      65             : This metric term demonstrates that a gridpoint distribution uniform in
      66             : \f$\xi\f$ will end up being compressed near \f$\xi=\pm1\f$. Suppose we
      67             : reparameterized the surface using the generalized logical coordinate
      68             : \f$\Xi\in[-1,1]\f$. We would find:
      69             : 
      70             : \f[g_{\xi,\xi}|_{\eta=0} = \frac{\Xi'^2}{(1+\Xi^2)^2}\f]
      71             : 
      72             : Ideally, we would like distances measured along a curvilinear surface to be
      73             : linear in the logical coordinates. We solve the differential equation and
      74             : obtain:
      75             : 
      76             : \f[\Xi = \tan(\xi\pi/4)\f]
      77             : 
      78             : These two parameterizations of the cubed sphere are known as the *equidistant*
      79             : and *equiangular* central projections of the cube onto the sphere. We now
      80             : summarize their usage in SpECTRE CoordinateMaps that have
      81             : `with_equiangular_map` as a specifiable parameter:
      82             : 
      83             : In the case where `with_equiangular_map` is `true`, we have the
      84             : equiangular coordinates
      85             : 
      86             : \f[\textrm{equiangular xi} : \Xi(\xi) = \textrm{tan}(\xi\pi/4)\f]
      87             : 
      88             : \f[\textrm{equiangular eta}  : \mathrm{H}(\eta) = \textrm{tan}(\eta\pi/4)\f]
      89             : 
      90             : with derivatives
      91             : 
      92             : \f[\Xi'(\xi) = \frac{\pi}{4}(1+\Xi^2)\f],
      93             : 
      94             : \f[\mathrm{H}'(\eta) = \frac{\pi}{4}(1+\mathrm{H}^2)\f]
      95             : 
      96             : In the case where `with_equiangular_map` is `false`, we have the equidistant
      97             : coordinates
      98             : 
      99             : \f[ \textrm{equidistant xi}  : \Xi = \xi\f]
     100             : 
     101             : \f[ \textrm{equidistant eta}  : \mathrm{H} = \eta\f]
     102             : 
     103             : with derivatives:
     104             : 
     105             : \f[\Xi'(\xi) = 1\f] \f[\mathrm{H}'(\eta) = 1\f]
     106             : 
     107             : ##  Projective Maps
     108             : 
     109             : The mapping for any convex quadrilateral can be obtained by bilinearly
     110             : interpolating between each vertex \f$\vec{x}_1, \vec{x}_2, \vec{x}_3\f$
     111             : and \f$\vec{x}_4\f$:
     112             : 
     113             : \f[\vec{x}(\xi,\eta) =
     114             : \frac{(1-\xi)(1-\eta)}{4}\vec{x}_1+
     115             : \frac{(1+\xi)(1-\eta)}{4}\vec{x}_2+
     116             : \frac{(1-\xi)(1+\eta)}{4}\vec{x}_3+
     117             : \frac{(1+\xi)(1+\eta)}{4}\vec{x}_4
     118             : \f]
     119             : 
     120             : In the case of a trapezoid where two of the sides are parallel, it is
     121             : appropriate to linearly interpolate along the parallel sides. However,
     122             : linearly interpolating between the two bases results in a less than
     123             : ideal gridpoint distribution. This happens in the case of SpECTRE's Frustum,
     124             : where the logical coordinate \f$\zeta\f$ interpolates between the bases.
     125             : 
     126             : \image html BilinearVProjective.png "Comparison of mappings. (Noah Veltman)"
     127             : 
     128             : As seen in Veltman's [Warp-Off]
     129             : (https://bl.ocks.org/veltman/8f5a157276b1dc18ce2fba1bc06dfb48), linear
     130             : interpolation between the two bases results in a uniformly spaced grid
     131             : between the bases of the frustum. This causes elements near the smaller base
     132             : to be longer in the direction normal to the base, and elements near the larger
     133             : base to be shorter in the direction normal to the base. We desire elements that
     134             : have roughly equal sizes along each of their dimensions.
     135             : 
     136             : We can redistribute the gridpoints in the \f$\zeta\f$ direction using a
     137             : projective map, moving more gridpoints toward the smaller base. We can also see
     138             : in the figure above that a projective map can be applied incorrectly, leaving
     139             : elements distorted at the opposite end. From this we can see that it is
     140             : important to control the degree of projection.
     141             : 
     142             : We adapt a technique from projective geometry to obtain the desired grid
     143             : spacing. The heart of the method lies in the fact that objects arranged in a
     144             : line at equal distances from one another will appear to converge as they
     145             : approach the horizon.
     146             : 
     147             : \image html ProjectionOntoPlane.png "Controlling the degree of projection."
     148             : 
     149             : The above diagram demonstrates how to obtain a nonlinearly parameterized
     150             : object (seen in red) from a linearly parameterized one (seen in purple).
     151             : This is done by lifting the linearly parameterized object into a higher
     152             : spatial dimension \f$w\f$, such that its projection onto the plane remains
     153             : unchanged. As seen above, \f$w_{\delta}\f$ controls the degree of projection
     154             : of one end of the object (purple) into a higher spatial dimension \f$w\f$.
     155             : In projective geometry, these points that exist in the higher dimension are
     156             : labeled with *homogeneous coordinates* \f$\tilde{x}, \tilde{y}, \tilde{z}, w\f$,
     157             : to distinguish them from the Cartesian coordinates that label points that exist
     158             : on the \f$w=1\f$ hyperplane, \f$x,y,z\f$. The resulting grid (seen in red) is
     159             : obtained by projecting back into the \f$w=1\f$ hyperplane. The Cartesian
     160             : coordinates are obtained by dividing each homogeneous coordinate of the
     161             : linearly parameterized object by its respective \f$w\f$ coordinate value.
     162             : 
     163             : The parametric equation for the purple object seen above in homogeneous
     164             : coordinates is:
     165             : \f[\begin{bmatrix}\tilde{x}\\\tilde{y}\\\tilde{z}\\w\\\end{bmatrix}=
     166             : \frac{1-\zeta}{2}\begin{bmatrix}x_1\\y_1\\z_1\\1\\\end{bmatrix}+
     167             : \frac{1+\zeta}{2}\begin{bmatrix}x_2w_{\delta}\\y_2w_{\delta}\\
     168             : z_2w_{\delta}\\w_{\delta}\\\end{bmatrix}\f]
     169             : 
     170             : The equation for the projected red object in Cartesian coordinates is
     171             : obtained by dividing by w:
     172             : \f[\vec{x}(\zeta) = \frac{1}{w(\zeta)}
     173             : \begin{bmatrix}
     174             : \tilde{x}(\zeta)\\
     175             : \tilde{y}(\zeta)\\
     176             : \tilde{z}(\zeta)\\
     177             : \end{bmatrix}\f]
     178             : 
     179             : We wish to cast our parametric equation for the surface into the form:
     180             : \f[\vec{x}(\zeta) =
     181             : \frac{1-\mathrm{Z}}{2}\vec{x}_1 + \frac{1+\mathrm{Z}}{2}\vec{x}_2\f]
     182             : for some appropriate choice of auxiliary variable `projective_zeta`
     183             : \f$ = \mathrm{Z}(\zeta)\f$. We would also like for \f$\mathrm{Z}\f$ to reduce to
     184             : \f$\zeta\f$ when \f$w_{\delta}\ = 1\f$.
     185             : 
     186             : Defining the auxiliary variables \f$w_{\pm} := w_{\delta}\pm 1\f$, the desired
     187             : \f$\mathrm{Z}(\zeta)\f$ is given by:
     188             : \f[\mathrm{Z} = \frac{w_- + \zeta w_+}
     189             : {w_+ + \zeta w_-}\f]
     190             : 
     191             : with derivative:
     192             : \f[\mathrm{Z}' = \frac{\partial\mathrm{Z}}{\partial \zeta} =
     193             : \frac{w_+^2 - w_-^2}{(w_+ + \zeta w_-)^2}\f]
     194             : 
     195             : For SpECTRE CoordinateMaps that have `projective_scale_factor` as a specifiable
     196             : parameter, the value \f$w_{\delta} = 1\f$ should be supplied in case the user
     197             : does not want to use projective scaling. If `auto_projective_scale_factor` is
     198             : set to `true`, the map will compute a value of \f$w_{\delta}\f$ that is
     199             : appropriate.