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.
|
