Miscellaneous improvements (#1230)

* Refactor primitives

* More refactoring

Author: juliohm

src/geometries/primitives/ball.jl

@@ -50,10 +50,10 @@ function (b::Ball{𝔼{2}})(ρ, φ)
 end
 
 function (b::Ball{𝔼{3}})(ρ, θ, φ)
-  T = numtype(lentype(b))
   if (ρ < 0 || ρ > 1) || (θ < 0 || θ > 1) || (φ < 0 || φ > 1)
     throw(DomainError((ρ, θ, φ), "b(ρ, θ, φ) is not defined for ρ, θ, φ outside [0, 1]³."))
   end
+  T = numtype(lentype(b))
   c = b.center
   r = b.radius
   ρ′ = T(ρ) * r

src/geometries/primitives/circle.jl

@@ -43,24 +43,13 @@ paramdim(::Type{<:Circle}) = 1
 
 plane(c::Circle) = c.plane
 
-center(c::Circle) = c.plane(0, 0)
-
 radius(c::Circle) = c.radius
 
-==(c₁::Circle, c₂::Circle) = c₁.plane == c₂.plane && c₁.radius == c₂.radius
+center(c::Circle) = plane(c)(0, 0)
+
+==(c₁::Circle, c₂::Circle) = plane(c₁) == plane(c₂) && radius(c₁) == radius(c₂)
 
 Base.isapprox(c₁::Circle, c₂::Circle; atol=atol(lentype(c₁)), kwargs...) =
-  isapprox(c₁.plane, c₂.plane; atol, kwargs...) && isapprox(c₁.radius, c₂.radius; atol, kwargs...)
+  isapprox(plane(c₁), plane(c₂); atol, kwargs...) && isapprox(radius(c₁), radius(c₂); atol, kwargs...)
 
-function (c::Circle)(φ)
-  if (φ < 0 || φ > 1)
-    throw(DomainError(φ, "c(φ) is not defined for φ outside [0, 1]."))
-  end
-  T = numtype(lentype(c))
-  r = c.radius
-  l = r
-  sφ, cφ = sincospi(2 * T(φ))
-  u = ustrip(l * cφ)
-  v = ustrip(l * sφ)
-  c.plane(u, v)
-end
+(c::Circle)(φ) = Disk(plane(c), radius(c))(1, φ)

src/geometries/primitives/cone.jl

@@ -26,9 +26,9 @@ base(c::Cone) = c.base
 
 apex(c::Cone) = c.apex
 
-height(c::Cone) = norm(center(base(c)) - apex(c))
+height(c::Cone) = height(boundary(c))
 
-halfangle(c::Cone) = atan(radius(base(c)), height(c))
+halfangle(c::Cone) = halfangle(boundary(c))
 
 ==(c₁::Cone, c₂::Cone) = boundary(c₁) == boundary(c₂)
 
@@ -39,7 +39,7 @@ function (c::Cone)(r, φ, h)
   if (r < 0 || r > 1) || (φ < 0 || φ > 1) || (h < 0 || h > 1)
     throw(DomainError((r, φ, h), "c(r, φ, h) is not defined for r, φ, h outside [0, 1]³."))
   end
-  a = c.apex
-  b = c.base(r, φ)
+  b = base(c)(r, φ)
+  a = apex(c)
   Segment(b, a)(h)
 end

src/geometries/primitives/conesurface.jl

@@ -26,17 +26,13 @@ base(c::ConeSurface) = c.base
 
 apex(c::ConeSurface) = c.apex
 
-==(c₁::ConeSurface, c₂::ConeSurface) = c₁.base == c₂.base && c₁.apex == c₂.apex
+height(c::ConeSurface) = norm(center(base(c)) - apex(c))
+
+halfangle(c::ConeSurface) = atan(radius(base(c)), height(c))
+
+==(c₁::ConeSurface, c₂::ConeSurface) = base(c₁) == base(c₂) && apex(c₁) == apex(c₂)
 
 Base.isapprox(c₁::ConeSurface, c₂::ConeSurface; atol=atol(lentype(c₁)), kwargs...) =
-  isapprox(c₁.base, c₂.base; atol, kwargs...) && isapprox(c₁.apex, c₂.apex; atol, kwargs...)
-
-function (c::ConeSurface)(φ, h)
-  if (φ < 0 || φ > 1) || (h < 0 || h > 1)
-    throw(DomainError((φ, h), "c(φ, h) is not defined for φ, h outside [0, 1]²."))
-  end
-  T = numtype(lentype(c))
-  a = c.apex
-  b = c.base(one(T), φ)
-  Segment(b, a)(h)
-end
+  isapprox(base(c₁), base(c₂); atol, kwargs...) && isapprox(apex(c₁), apex(c₂); atol, kwargs...)
+
+(c::ConeSurface)(φ, h) = Cone(base(c), apex(c))(1, φ, h)

src/geometries/primitives/cylinder.jl

@@ -62,6 +62,17 @@ top(c::Cylinder) = c.top
 
 radius(c::Cylinder) = c.radius
 
+axis(c::Cylinder) = axis(boundary(c))
+
+isright(c::Cylinder) = isright(boundary(c))
+
+hasintersectingplanes(c::Cylinder) = hasintersectingplanes(boundary(c))
+
+==(c₁::Cylinder, c₂::Cylinder) = boundary(c₁) == boundary(c₂)
+
+Base.isapprox(c₁::Cylinder, c₂::Cylinder; atol=atol(lentype(c₁)), kwargs...) =
+  isapprox(boundary(c₁), boundary(c₂); atol, kwargs...)
+
 function (c::Cylinder)(ρ, φ, z)
   if (ρ < 0 || ρ > 1) || (φ < 0 || φ > 1) || (z < 0 || z > 1)
     throw(DomainError((ρ, φ, z), "c(ρ, φ, z) is not defined for ρ, φ, z outside [0, 1]³."))
@@ -89,18 +100,3 @@ function (c::Cylinder)(ρ, φ, z)
   s = Segment(l ∩ b, l ∩ t)
   s(T(z))
 end
-
-# ----------------------------------------------
-# forward methods to boundary (CylinderSurface)
-# ----------------------------------------------
-
-axis(c::Cylinder) = axis(boundary(c))
-
-isright(c::Cylinder) = isright(boundary(c))
-
-hasintersectingplanes(c::Cylinder) = hasintersectingplanes(boundary(c))
-
-==(c₁::Cylinder, c₂::Cylinder) = boundary(c₁) == boundary(c₂)
-
-Base.isapprox(c₁::Cylinder, c₂::Cylinder; atol=atol(lentype(c₁)), kwargs...) =
-  isapprox(boundary(c₁), boundary(c₂); atol, kwargs...)

src/geometries/primitives/disk.jl

@@ -22,22 +22,22 @@ paramdim(::Type{<:Disk}) = 2
 
 plane(d::Disk) = d.plane
 
-center(d::Disk) = d.plane(0, 0)
-
 radius(d::Disk) = d.radius
 
-normal(d::Disk) = normal(d.plane)
+center(d::Disk) = center(boundary(d))
+
+normal(d::Disk) = normal(plane(d))
 
-==(d₁::Disk, d₂::Disk) = d₁.plane == d₂.plane && d₁.radius == d₂.radius
+==(d₁::Disk, d₂::Disk) = boundary(d₁) == boundary(d₂)
 
 Base.isapprox(d₁::Disk, d₂::Disk; atol=atol(lentype(d₁)), kwargs...) =
-  isapprox(d₁.plane, d₂.plane; atol, kwargs...) && isapprox(d₁.radius, d₂.radius; atol, kwargs...)
+  isapprox(boundary(d₁), boundary(d₂); atol, kwargs...)
 
 function (d::Disk)(ρ, φ)
-  T = numtype(lentype(d))
   if (ρ < 0 || ρ > 1) || (φ < 0 || φ > 1)
     throw(DomainError((ρ, φ), "d(ρ, φ) is not defined for ρ, φ outside [0, 1]²."))
   end
+  T = numtype(lentype(d))
   r = d.radius
   l = T(ρ) * r
   sφ, cφ = sincospi(2 * T(φ))

src/geometries/primitives/ellipsoid.jl

@@ -42,10 +42,10 @@ function Base.isapprox(e₁::Ellipsoid, e₂::Ellipsoid; atol=atol(lentype(e₁)
 end
 
 function (e::Ellipsoid)(θ, φ)
-  T = numtype(lentype(e))
   if (θ < 0 || θ > 1) || (φ < 0 || φ > 1)
     throw(DomainError((θ, φ), "e(θ, φ) is not defined for θ, φ outside [0, 1]²."))
   end
+  T = numtype(lentype(e))
   r = e.radii
   c = e.center
   R = e.rotation

src/geometries/primitives/frustumsurface.jl

@@ -36,35 +36,9 @@ height(f::FrustumSurface) = norm(center(bottom(f)) - center(top(f)))
 
 axis(f::FrustumSurface) = Line(center(bottom(f)), center(top(f)))
 
-==(f₁::FrustumSurface, f₂::FrustumSurface) = f₁.bot == f₂.bot && f₁.top == f₂.top
+==(f₁::FrustumSurface, f₂::FrustumSurface) = bottom(f₁) == bottom(f₂) && top(f₁) == top(f₂)
 
 Base.isapprox(f₁::FrustumSurface, f₂::FrustumSurface; atol=atol(lentype(f₁)), kwargs...) =
-  isapprox(f₁.bot, f₂.bot; atol, kwargs...) && isapprox(f₁.top, f₂.top; atol, kwargs...)
+  isapprox(bottom(f₁), bottom(f₂); atol, kwargs...) && isapprox(top(f₁), top(f₂); atol, kwargs...)
 
-function (f::FrustumSurface)(φ, z)
-  ℒ = lentype(f)
-  T = numtype(ℒ)
-  if (φ < 0 || φ > 1) || (z < 0 || z > 1)
-    throw(DomainError((φ, z), "f(φ, z) is not defined for φ, z outside [0, 1]²."))
-  end
-  rb = radius(bottom(f))
-  rt = radius(top(f))
-  a = axis(f)
-  d = a(1) - a(0)
-  l = norm(d)
-
-  # rotation to align z axis with cylinder axis
-  Q = urotbetween(d, Vec(zero(ℒ), zero(ℒ), oneunit(ℒ)))
-
-  # scale coordinates
-  φₛ = 2T(π) * φ
-  zₛ = z * l
-
-  # local coordinates, that will be transformed with rotation and position of the FrustumSurface
-  x = cos(φₛ) * (rb * (l - zₛ) + rt * zₛ) / l
-  y = sin(φₛ) * (rb * (l - zₛ) + rt * zₛ) / l
-  z = zₛ
-  p = Vec(x, y, z)
-
-  center(bottom(f)) + Q' * p
-end
+(f::FrustumSurface)(φ, z) = Frustum(bottom(f), top(f))(1, φ, z)