Add inverse cube sphere warping

Author: lcw

src/gridgenerators.jl

@@ -278,6 +278,26 @@ function cubespherewarp(point)
     return point
 end
 
+function cubesphereunwarp(point)
+    # Put the points in reverse magnitude order
+    p = sortperm(abs.(point))
+    point = point[p]
+
+    # Convert to angles
+    ξ = 4atan(point[2]/point[3])/π
+    η = 4atan(point[1]/point[3])/π
+    R = sign(point[3])*hypot(point...)
+
+    x = R
+    y = R * ξ
+    z = R * η
+
+    # Compute the new points and unpermute
+    point = SVector(z, y, x)[sortperm(p)]
+
+    return point
+end
+
 function cubedshellconnectivity(referencecell::LobattoCell,
                                 R::Real,
                                 ncells::Integer;

test/gridgenerators.jl

@@ -239,4 +239,18 @@ end
             end
         end
     end
+
+    for p in (
+            SVector(+0.4, -1.00, +0.3),
+            SVector(+0.5, +1.00, -0.2),
+            SVector(-0.5, +1.00, +0.2),
+            SVector(+1.0, -0.70, -0.2),
+            SVector(-1.0, -0.70, +0.2),
+            SVector(+0.3, -0.75, -1.0),
+            SVector(-0.3, -0.75, +1.0),
+        )
+        r = Bennu.cubesphereunwarp(p)
+        p̂ = Bennu.cubespherewarp(r)
+        @test p ≈ p̂
+    end
 end