Update parameterization implementation of Cylinder and `CylinderSur…

…face` & Simplification of `Sphere` parameterization (#1098) * Update parameterization implementation of 'Cylinder' * Fix typo * Update parameterization implementation of 'CylinderSurface' * Update parameterization implementation of 'Sphere' * Apply suggestions * Apply suggestions * Add more tests * Apply suggestions * Update code * Apply suggestions
JuliaGeometry · Oct 4, 2024 · ca2b918 · ca2b918
1 parent 77b2384
commit ca2b918
Show file tree

Hide file tree

Showing 4 changed files with 32 additions and 62 deletions.
diff --git a/src/geometries/primitives/cylinder.jl b/src/geometries/primitives/cylinder.jl
@@ -83,16 +83,17 @@ function (c::Cylinder)(ρ, φ, z)
   a = axis(c)
   d = a(T(1)) - a(T(0))
   h = norm(d)
-  o = b(0, 0)
+  o = b(T(0), T(0))
 
   # rotation to align z axis with cylinder axis
   Q = urotbetween(Vec(zero(ℒ), zero(ℒ), oneunit(ℒ)), d)
 
   # project a parametric segment between the top and bottom planes
-  lsφ, lcφ = T(ρ) * r .* sincospi(2 * T(φ))
-  p₁ = o + Q * Vec(lcφ, lsφ, zero(ℒ))
-  p₂ = o + Q * Vec(lcφ, lsφ, h)
-  l = Line(p₁, p₂)
+  ρ′ = T(ρ) * r
+  φ′ = T(φ) * 2 * T(π) * u"rad"
+  p₁ = Point(convert(crs(c), Cylindrical(ρ′, φ′, zero(ℒ))))
+  p₂ = Point(convert(crs(c), Cylindrical(ρ′, φ′, h)))
+  l = Line(p₁, p₂) |> Affine(Q, to(o))
   s = Segment(l ∩ b, l ∩ t)
   s(T(z))
 end
diff --git a/src/geometries/primitives/cylindersurface.jl b/src/geometries/primitives/cylindersurface.jl
@@ -83,39 +83,7 @@ Base.isapprox(c₁::CylinderSurface, c₂::CylinderSurface; atol=atol(lentype(c
   isapprox(c₁.top, c₂.top; atol, kwargs...) &&
   isapprox(c₁.radius, c₂.radius; atol, kwargs...)
 
-function (c::CylinderSurface)(φ, z)
-  ℒ = lentype(c)
-  T = numtype(ℒ)
-  if (φ < 0 || φ > 1) || (z < 0 || z > 1)
-    throw(DomainError((φ, z), "c(φ, z) is not defined for φ, z outside [0, 1]²."))
-  end
-  t = top(c)
-  b = bottom(c)
-  r = radius(c)
-  a = axis(c)
-  d = a(T(1)) - a(T(0))
-  h = norm(d)
-  o = centroid(c)
-
-  # rotation to align z axis with cylinder axis
-  Q = urotbetween(d, Vec(zero(ℒ), zero(ℒ), oneunit(ℒ)))
-
-  # new normals of planes in new rotated system
-  nᵦ = Q * normal(b)
-  nₜ = Q * normal(t)
-
-  # given cylindrical coordinates (r*cos(φ), r*sin(φ), z) and the
-  # equation of the plane, we can solve for z and find all points
-  # along the ellipse obtained by intersection
-  rsφ, rcφ = r .* sincospi(2 * T(φ))
-  zᵦ = -h / 2 - (rcφ * nᵦ[1] + rsφ * nᵦ[2]) / nᵦ[3]
-  zₜ = +h / 2 - (rcφ * nₜ[1] + rsφ * nₜ[2]) / nₜ[3]
-  pᵦ = Point(rcφ, rsφ, zᵦ)
-  pₜ = Point(rcφ, rsφ, zₜ)
-
-  p = pᵦ + T(z) * (pₜ - pᵦ)
-  o + Q' * to(p)
-end
+(c::CylinderSurface)(φ, z) = Cylinder(bottom(c), top(c), radius(c))(1, φ, z)
 
 function hasintersectingplanes(c::CylinderSurface)
   x = c.bot ∩ c.top

diff --git a/src/geometries/primitives/sphere.jl b/src/geometries/primitives/sphere.jl
@@ -73,27 +73,4 @@ radius(s::Sphere) = s.radius
 Base.isapprox(s₁::Sphere, s₂::Sphere; atol=atol(lentype(s₁)), kwargs...) =
   isapprox(s₁.center, s₂.center; atol, kwargs...) && isapprox(s₁.radius, s₂.radius; atol, kwargs...)
 
-function (s::Sphere{𝔼{2}})(φ)
-  T = numtype(lentype(s))
-  if (φ < 0 || φ > 1)
-    throw(DomainError(φ, "s(φ) is not defined for φ outside [0, 1]."))
-  end
-  c = s.center
-  r = s.radius
-  φ′ = T(φ) * 2 * T(π) * u"rad"
-  p = Point(convert(crs(s), Polar(r, φ′)))
-  p + to(c)
-end
-
-function (s::Sphere{𝔼{3}})(θ, φ)
-  T = numtype(lentype(s))
-  if (θ < 0 || θ > 1) || (φ < 0 || φ > 1)
-    throw(DomainError((θ, φ), "s(θ, φ) is not defined for θ, φ outside [0, 1]²."))
-  end
-  c = s.center
-  r = s.radius
-  θ′ = T(θ) * T(π) * u"rad"
-  φ′ = T(φ) * 2 * T(π) * u"rad"
-  p = Point(convert(crs(s), Spherical(r, θ′, φ′)))
-  p + to(c)
-end
+(s::Sphere)(uv...) = Ball(center(s), radius(s))(1, uv...)
diff --git a/test/primitives.jl b/test/primitives.jl
@@ -594,6 +594,12 @@ end
   @test b(T(0), T(0)) ≈ cart(0, 0)
   @test b(T(1), T(0)) ≈ cart(2, 0)
 
+  # machine type is preserved in parameterization
+  b = Ball(cart(0, 0), T(2))
+  @test Meshes.lentype(b(0, 0)) == ℳ
+  @test Meshes.lentype(b(0.0, 0.0)) == ℳ
+  @test Meshes.lentype(b(0.0f0, 0.0f0)) == ℳ
+
   b = Ball(cart(7, 7), T(1.5))
   ps = b.(1, rand(T, 100))
   all(∈(b), ps)
@@ -709,6 +715,12 @@ end
   @test s(T(0), T(0)) ≈ cart(0, 0, 2)
   @test s(T(0.5), T(0.5)) ≈ cart(-2, 0, 0)
 
+  # machine type is preserved in parameterization
+  s = Sphere(cart(0, 0), T(2))
+  @test Meshes.lentype(s(0)) == ℳ
+  @test Meshes.lentype(s(0.0)) == ℳ
+  @test Meshes.lentype(s(0.0f0)) == ℳ
+
   s = Sphere(cart(0, 0, 0), T(1))
   @test sprint(show, s) == "Sphere(center: (x: 0.0 m, y: 0.0 m, z: 0.0 m), radius: 1.0 m)"
   if T === Float32
@@ -918,6 +930,12 @@ end
   @test cart(0, 0, 1.001) ∉ c
   @test cart(1, 1, 1) ∉ c
 
+  # machine type is preserved in parameterization
+  c = Cylinder(T(1))
+  @test Meshes.lentype(c(0, 0, 0)) == ℳ
+  @test Meshes.lentype(c(0.0, 0.0, 0.0)) == ℳ
+  @test Meshes.lentype(c(0.0f0, 0.0f0, 0.0f0)) == ℳ
+
   c = Cylinder(cart(0, 0, 0), cart(0, 0, 1), T(1))
   @test sprint(show, c) ==
         "Cylinder(bot: Plane(p: (x: 0.0 m, y: 0.0 m, z: 0.0 m), u: (1.0 m, -0.0 m, -0.0 m), v: (-0.0 m, 1.0 m, -0.0 m)), top: Plane(p: (x: 0.0 m, y: 0.0 m, z: 1.0 m), u: (1.0 m, -0.0 m, -0.0 m), v: (-0.0 m, 1.0 m, -0.0 m)), radius: 1.0 m)"
@@ -974,6 +992,12 @@ end
   @test c(T(0), T(1)) ≈ cart(1, 0, 1)
   @test c(T(0.5), T(1)) ≈ cart(-1, 0, 1)
 
+  # machine type is preserved in parameterization
+  c = CylinderSurface(T(1))
+  @test Meshes.lentype(c(0, 0)) == ℳ
+  @test Meshes.lentype(c(0.0, 0.0)) == ℳ
+  @test Meshes.lentype(c(0.0f0, 0.0f0)) == ℳ
+
   c = CylinderSurface(1.0)
   @test Meshes.lentype(c) == Meshes.Met{Float64}
   c = CylinderSurface(1.0f0)