Miscellaneous improvements

src/discretization/regular.jl

@@ -27,6 +27,12 @@ function discretize(geometry::Geometry, method::RegularDiscretization)
   end
 end
 
+# box is trivial to discretize into a regular grid
+function discretize(box::Box, method::RegularDiscretization)
+  sz = fitdims(method.sizes, paramdim(box))
+  RegularGrid(extrema(box)..., dims=sz)
+end
+
 # triangle is parametrized with barycentric coordinates, so we can't rely on regular sampling
 discretize(triangle::Triangle, method::RegularDiscretization) = discretizewithinbox(triangle, method)
 
@@ -180,12 +186,3 @@ function _appendpoles(tg, d, ccw)
 
   SimpleTopology([trunk; north; south])
 end
-
-# --------------
-# SPECIAL CASES
-# --------------
-
-function discretize(box::Box, method::RegularDiscretization)
-  sz = fitdims(method.sizes, paramdim(box))
-  RegularGrid(extrema(box)..., dims=sz)
-end

src/geometries/primitives/ball.jl

@@ -37,10 +37,10 @@ Base.isapprox(b₁::Ball, b₂::Ball; atol=atol(lentype(b₁)), kwargs...) =
   isapprox(b₁.center, b₂.center; atol, kwargs...) && isapprox(b₁.radius, b₂.radius; atol, kwargs...)
 
 function (b::Ball{𝔼{2}})(ρ, φ)
-  T = numtype(lentype(b))
   if (ρ < 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/bezier.jl

@@ -7,12 +7,12 @@
 
 A recursive Bézier curve with control points `points`.
 See <https://en.wikipedia.org/wiki/Bézier_curve>.
-A point on the curve `b` can be evaluated by calling
-`b(t)` with `t` between `0` and `1`.
-The evaluation method defaults to DeCasteljau's algorithm
-for accurate evaluation. Horner's method, faster with a
-large number of points but less precise, can be used via
-`b(t, Horner())`.
+
+A point on the Bézier curve `b` can be evaluated by calling
+`b(t)` with `t` between `0` and `1`. The evaluation method
+defaults to DeCasteljau's algorithm for accurate evaluation.
+Horner's method, faster with a large number of points but
+less precise, can be used via `b(t, Horner())`.
 
 ## Examples
 
@@ -86,10 +86,10 @@ end
 # Horner's rule recursively reconstructs B from a sequence bᵢ
 # with bₙ = aₙ and bᵢ₋₁ = aᵢ₋₁ + bᵢ * t until b₀ = B.
 function (curve::BezierCurve)(t, ::Horner)
-  T = numtype(lentype(curve))
   if t < 0 || t > 1
     throw(DomainError(t, "b(t) is not defined for t outside [0, 1]."))
   end
+  T = numtype(lentype(curve))
   cs = curve.controls
   t̄ = one(T) - t
   n = degree(curve)

src/geometries/primitives/circle.jl

@@ -53,10 +53,10 @@ 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...)
 
 function (c::Circle)(φ)
-  T = numtype(lentype(c))
   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(φ))

src/geometries/primitives/conesurface.jl

@@ -32,10 +32,10 @@ Base.isapprox(c₁::ConeSurface, c₂::ConeSurface; atol=atol(lentype(c₁)), kw
   isapprox(c₁.base, c₂.base; atol, kwargs...) && isapprox(c₁.apex, c₂.apex; atol, kwargs...)
 
 function (c::ConeSurface)(φ, h)
-  T = numtype(lentype(c))
   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)

src/geometries/primitives/cylinder.jl

@@ -22,7 +22,9 @@ with `start` and `finish` end points.
 
 Finally, construct a right vertical circular cylinder with given `radius`.
 
-See <https://en.wikipedia.org/wiki/Cylinder>. 
+See <https://en.wikipedia.org/wiki/Cylinder>.
+
+See also [`CylinderSurface`](@ref).
 """
 struct Cylinder{C<:CRS,P<:Plane{C},ℒ<:Len} <: Primitive{𝔼{3},C}
   bot::P
@@ -54,46 +56,51 @@ end
 
 paramdim(::Type{<:Cylinder}) = 3
 
-radius(c::Cylinder) = c.radius
-
 bottom(c::Cylinder) = c.bot
 
 top(c::Cylinder) = c.top
 
-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...)
+radius(c::Cylinder) = c.radius
 
 function (c::Cylinder)(ρ, φ, z)
-  ℒ = lentype(c)
-  T = numtype(ℒ)
   if (ρ < 0 || ρ > 1) || (φ < 0 || φ > 1) || (z < 0 || z > 1)
     throw(DomainError((ρ, φ, z), "c(ρ, φ, z) is not defined for ρ, φ, z outside [0, 1]³."))
   end
+  ℒ = lentype(c)
+  T = numtype(ℒ)
   t = top(c)
   b = bottom(c)
   r = radius(c)
   a = axis(c)
   d = a(T(1)) - a(T(0))
-  h = norm(d)
-  o = b(T(0), T(0))
 
   # rotation to align z axis with cylinder axis
-  Q = urotbetween(Vec(zero(ℒ), zero(ℒ), oneunit(ℒ)), d)
+  R = urotbetween(Vec(zero(ℒ), zero(ℒ), oneunit(ℒ)), d)
+
+  # offset to translate cylinder to final position
+  o = to(b(T(0), T(0)))
 
   # project a parametric segment between the top and bottom planes
   ρ′ = 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))
+  p₂ = Point(convert(crs(c), Cylindrical(ρ′, φ′, norm(d))))
+  l = Line(p₁, p₂) |> Affine(R, o)
   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/cylindersurface.jl

@@ -22,7 +22,9 @@ with `start` and `finish` end points.
 
 Finally, construct a right vertical circular cylinder surface with given `radius`.
 
-See <https://en.wikipedia.org/wiki/Cylinder>. 
+See <https://en.wikipedia.org/wiki/Cylinder>.
+
+See also [`Cylinder`](@ref).
 """
 struct CylinderSurface{C<:CRS,P<:Plane{C},ℒ<:Len} <: Primitive{𝔼{3},C}
   bot::P
@@ -54,38 +56,35 @@ end
 
 paramdim(::Type{<:CylinderSurface}) = 2
 
-radius(c::CylinderSurface) = c.radius
-
 bottom(c::CylinderSurface) = c.bot
 
 top(c::CylinderSurface) = c.top
 
-axis(c::CylinderSurface) = Line(c.bot(0, 0), c.top(0, 0))
+radius(c::CylinderSurface) = c.radius
+
+axis(c::CylinderSurface) = Line(bottom(c)(0, 0), top(c)(0, 0))
 
+# cylinder is right if axis is aligned with plane normals
 function isright(c::CylinderSurface)
-  ℒ = lentype(c)
-  T = numtype(ℒ)
-  # cylinder is right if axis
-  # is aligned with plane normals
+  T = numtype(lentype(c))
   a = axis(c)
   d = a(T(1)) - a(T(0))
-  v = normal(c.bot)
-  w = normal(c.top)
-  isparallelv = isapproxzero(norm(d × v))
-  isparallelw = isapproxzero(norm(d × w))
-  isparallelv && isparallelw
+  u = normal(bottom(c))
+  v = normal(top(c))
+  isapproxzero(norm(d × u)) && isapproxzero(norm(d × v))
 end
 
-==(c₁::CylinderSurface, c₂::CylinderSurface) = c₁.bot == c₂.bot && c₁.top == c₂.top && c₁.radius == c₂.radius
+function hasintersectingplanes(c::CylinderSurface)
+  l = bottom(c) ∩ top(c)
+  !isnothing(l) && evaluate(Euclidean(), axis(c), l) < radius(c)
+end
+
+==(c₁::CylinderSurface, c₂::CylinderSurface) =
+  bottom(c₁) == bottom(c₂) && top(c₁) == top(c₂) && radius(c₁) == radius(c₂)
 
 Base.isapprox(c₁::CylinderSurface, c₂::CylinderSurface; atol=atol(lentype(c₁)), kwargs...) =
-  isapprox(c₁.bot, c₂.bot; atol, kwargs...) &&
-  isapprox(c₁.top, c₂.top; atol, kwargs...) &&
-  isapprox(c₁.radius, c₂.radius; atol, kwargs...)
+  isapprox(bottom(c₁), bottom(c₂); atol, kwargs...) &&
+  isapprox(top(c₁), top(c₂); atol, kwargs...) &&
+  isapprox(radius(c₁), radius(c₂); atol, kwargs...)
 
 (c::CylinderSurface)(φ, z) = Cylinder(bottom(c), top(c), radius(c))(1, φ, z)
-
-function hasintersectingplanes(c::CylinderSurface)
-  x = c.bot ∩ c.top
-  !isnothing(x) && evaluate(Euclidean(), axis(c), x) < c.radius
-end

src/pointification.jl

@@ -37,6 +37,8 @@
 
 pointify(c::CylinderSurface) = _rsample(c)
 
+pointify(c::ConeSurface) = _rsample(c)
+
 pointify(p::PolyArea) = vertices(p)
 
 pointify(r::Ring) = vertices(r)