Bottom top fixes (#1332)

* Fixes to botom and top of Cylinder

* Fix Chain parametrization

Author: juliohm

src/boundary.jl

@@ -122,7 +122,7 @@ boundary(::Circle) = nothing
 
 embedboundary(c::Circle) = c
 
-boundary(c::Cylinder) = CylinderSurface(bottom(c), top(c), radius(c))
+boundary(c::Cylinder) = CylinderSurface(plane(bottom(c)), plane(top(c)), radius(c))
 
 embedboundary(c::Cylinder) = boundary(c)
 

src/geometries/polytopes.jl

@@ -157,20 +157,22 @@ function angles(c::Chain)
 end
 
 function (c::Chain)(t)
-  if t < 0 || t > 1
-    throw(DomainError(t, "c(t) is not defined for t outside [0, 1]."))
-  end
   segs = segments(c)
   sums = cumsum(measure.(segs))
   sums /= last(sums)
-  # find k such that sums[k] ≤ t < sums[k + 1]
-  k = searchsortedfirst(sums, t) - 1
-  # select segment s at index k
-  s, _ = iterate(segs, k)
+  # find first k such that t ≤ sums[k]
+  k = searchsortedfirst(sums, clamp(t, zero(t), one(t)))
   # reparametrization of t within s
-  ∑ₖ = iszero(k) ? zero(eltype(sums)) : sums[k]
-  ∑ₖ₊₁ = sums[k + 1]
-  s((t - ∑ₖ) / (∑ₖ₊₁ - ∑ₖ))
+  if isone(k)
+    s = first(segs)
+    Σₖ = sums[k]
+    s(t / Σₖ)
+  else
+    s = first(Iterators.drop(segs, k - 1))
+    Σₖ = sums[k]
+    Σₖ₋₁ = sums[k - 1]
+    s((t - Σₖ₋₁) / (Σₖ - Σₖ₋₁))
+  end
 end
 
 # implementations of Chain

src/geometries/primitives/cylinder.jl

@@ -56,9 +56,13 @@ end
 
 paramdim(::Type{<:Cylinder}) = 3
 
-bottom(c::Cylinder) = c.bot
+bottom(c::Cylinder) = Disk(c.bot, bottomradius(c))
 
-top(c::Cylinder) = c.top
+top(c::Cylinder) = Disk(c.top, topradius(c))
+
+bottomradius(c::Cylinder) = norm(c(1, 0, 0) - c(0, 0, 0))
+
+topradius(c::Cylinder) = norm(c(1, 0, 1) - c(0, 0, 1))
 
 radius(c::Cylinder) = c.radius
 
@@ -76,10 +80,12 @@ Base.isapprox(c₁::Cylinder, c₂::Cylinder; atol=atol(lentype(c₁)), kwargs..
 function (c::Cylinder)(ρ, φ, z)
   ℒ = lentype(c)
   T = numtype(ℒ)
-  t = top(c)
-  b = bottom(c)
-  r = radius(c)
-  a = axis(c)
+  C = crs(c)
+  D = datum(C)
+  b = c.bot
+  t = c.top
+  r = c.radius
+  a = Line(b(0, 0), t(0, 0))
   d = a(T(1)) - a(T(0))
 
   # rotation to align z axis with cylinder axis
@@ -91,8 +97,8 @@ function (c::Cylinder)(ρ, φ, z)
   # 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(ρ′, φ′, norm(d))))
+  p₁ = Point(convert(C, Cylindrical{D}(ρ′, φ′, zero(ℒ))))
+  p₂ = Point(convert(C, Cylindrical{D}(ρ′, φ′, norm(d))))
   l = Line(p₁, p₂) |> Affine(R, o)
   s = Segment(l ∩ b, l ∩ t)
   s(T(z))

src/geometries/primitives/cylindersurface.jl

@@ -56,35 +56,38 @@ end
 
 paramdim(::Type{<:CylinderSurface}) = 2
 
-bottom(c::CylinderSurface) = c.bot
+bottom(c::CylinderSurface) = Disk(c.bot, bottomradius(c))
 
-top(c::CylinderSurface) = c.top
+top(c::CylinderSurface) = Disk(c.top, topradius(c))
+
+bottomradius(c::CylinderSurface) = bottomradius(Cylinder(c.bot, c.top, c.radius))
+
+topradius(c::CylinderSurface) = topradius(Cylinder(c.bot, c.top, c.radius))
 
 radius(c::CylinderSurface) = c.radius
 
-axis(c::CylinderSurface) = Line(bottom(c)(0, 0), top(c)(0, 0))
+axis(c::CylinderSurface) = Line(c.bot(0, 0), c.top(0, 0))
 
 # cylinder is right if axis is aligned with plane normals
 function isright(c::CylinderSurface)
   T = numtype(lentype(c))
   a = axis(c)
   d = a(T(1)) - a(T(0))
-  u = normal(bottom(c))
-  v = normal(top(c))
+  u = normal(c.bot)
+  v = normal(c.top)
   isapproxzero(norm(d × u)) && isapproxzero(norm(d × v))
 end
 
 function hasintersectingplanes(c::CylinderSurface)
-  l = bottom(c) ∩ top(c)
-  !isnothing(l) && evaluate(Euclidean(), axis(c), l) < radius(c)
+  l = c.bot ∩ c.top
+  !isnothing(l) && evaluate(Euclidean(), axis(c), l) < c.radius
 end
 
-==(c₁::CylinderSurface, c₂::CylinderSurface) =
-  bottom(c₁) == bottom(c₂) && top(c₁) == top(c₂) && radius(c₁) == radius(c₂)
+==(c₁::CylinderSurface, c₂::CylinderSurface) = c₁.bot == c₂.bot && c₁.top == c₂.top && c₁.radius == c₂.radius
 
 Base.isapprox(c₁::CylinderSurface, c₂::CylinderSurface; atol=atol(lentype(c₁)), kwargs...) =
-  isapprox(bottom(c₁), bottom(c₂); atol, kwargs...) &&
-  isapprox(top(c₁), top(c₂); atol, kwargs...) &&
-  isapprox(radius(c₁), radius(c₂); atol, kwargs...)
+  isapprox(c₁.bot, c₂.bot; atol, kwargs...) &&
+  isapprox(c₁.top, c₂.top; atol, kwargs...) &&
+  isapprox(c₁.radius, c₂.radius; atol, kwargs...)
 
-(c::CylinderSurface)(φ, z) = Cylinder(bottom(c), top(c), radius(c))(1, φ, z)
+(c::CylinderSurface)(φ, z) = Cylinder(c.bot, c.top, c.radius)(1, φ, z)

test/primitives.jl

@@ -863,8 +863,8 @@ end
   @test crs(c) <: Cartesian{NoDatum}
   @test Meshes.lentype(c) == ℳ
   @test radius(c) == T(5) * u"m"
-  @test bottom(c) == Plane(cart(1, 2, 3), vector(0, 0, 1))
-  @test top(c) == Plane(cart(4, 5, 6), vector(0, 0, 1))
+  @test bottom(c) == Disk(Plane(cart(1, 2, 3), vector(0, 0, 1)), Meshes.bottomradius(c))
+  @test top(c) == Disk(Plane(cart(4, 5, 6), vector(0, 0, 1)), Meshes.topradius(c))
   @test axis(c) == Line(cart(1, 2, 3), cart(4, 5, 6))
   @test !isright(c)
   @test measure(c) == volume(c) ≈ T(5)^2 * pi * T(3) * sqrt(T(3)) * u"m^3"
@@ -897,8 +897,8 @@ end
 
   c = Cylinder(cart(0, 0, 0), cart(0, 0, 1), T(1))
   @test radius(c) == T(1) * u"m"
-  @test bottom(c) == Plane(cart(0, 0, 0), vector(0, 0, 1))
-  @test top(c) == Plane(cart(0, 0, 1), vector(0, 0, 1))
+  @test bottom(c) == Disk(Plane(cart(0, 0, 0), vector(0, 0, 1)), T(1))
+  @test top(c) == Disk(Plane(cart(0, 0, 1), vector(0, 0, 1)), T(1))
   @test centroid(c) == cart(0.0, 0.0, 0.5)
   @test axis(c) == Line(cart(0, 0, 0), cart(0, 0, 1))
   @test isright(c)
@@ -943,8 +943,8 @@ end
   @test crs(c) <: Cartesian{NoDatum}
   @test Meshes.lentype(c) == ℳ
   @test radius(c) == T(2) * u"m"
-  @test bottom(c) == Plane(cart(0, 0, 0), vector(0, 0, 1))
-  @test top(c) == Plane(cart(0, 0, 1), vector(0, 0, 1))
+  @test bottom(c) == Disk(Plane(cart(0, 0, 0), vector(0, 0, 1)), T(2))
+  @test top(c) == Disk(Plane(cart(0, 0, 1), vector(0, 0, 1)), T(2))
   @test centroid(c) == cart(0.0, 0.0, 0.5)
   @test axis(c) == Line(cart(0, 0, 0), cart(0, 0, 1))
   @test isright(c)