From 2f9d0d0212a37b5868bd81cc93cf56e4ede6a07e Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 12:24:08 -0500
Subject: [PATCH 01/15] Use hcubature with buffer for Rope

---
 src/specializations/Rope.jl | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index e917708b..0b4f27a5 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -35,3 +35,19 @@ function integral(
     # Convert the Rope into Segments, sum the integrals of those
     return sum(segment -> integral(f, segment, rule; kwargs...), Meshes.segments(rope))
 end
+
+function integral(
+    f,
+    rope::Meshes.Rope,
+    rule::HAdaptiveCubature;
+    FP::Type{T} = Float64,
+    kwargs...
+)
+    # Append a buffer to the given rule
+    buffer = HCubature.hcubature_buffer(f, _zeros(FP, 1), _ones(FP, 2))
+    rule = HAdaptiveCubature(rule.kwargs..., buffer = buffer)
+
+    # Convert the Rope into Segments, sum the integrals of those
+    _subintegral(seg) = _integral(f, seg, rule; FP = FP, rule.kwargs...)
+    return sum(_subintegral, Meshes.segments(rope))
+end

From c0bff58d4838960d6a6fa7a688b5ea8ca574c707 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 12:54:47 -0500
Subject: [PATCH 02/15] Add type param

---
 src/specializations/Rope.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index 0b4f27a5..064490b5 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -42,12 +42,12 @@ function integral(
     rule::HAdaptiveCubature;
     FP::Type{T} = Float64,
     kwargs...
-)
+) where {T <: AbstractFloat}
     # Append a buffer to the given rule
     buffer = HCubature.hcubature_buffer(f, _zeros(FP, 1), _ones(FP, 2))
     rule = HAdaptiveCubature(rule.kwargs..., buffer = buffer)
 
     # Convert the Rope into Segments, sum the integrals of those
-    _subintegral(seg) = _integral(f, seg, rule; FP = FP, rule.kwargs...)
+    _subintegral(seg) = _integral(f, seg, rule; FP = FP, kwargs...)
     return sum(_subintegral, Meshes.segments(rope))
 end

From 179da14f3cabc814a12fc12769c091b0f9ef86f9 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 13:00:41 -0500
Subject: [PATCH 03/15] Fix typo

---
 src/specializations/Rope.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index 064490b5..594f1a6a 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -44,7 +44,7 @@ function integral(
     kwargs...
 ) where {T <: AbstractFloat}
     # Append a buffer to the given rule
-    buffer = HCubature.hcubature_buffer(f, _zeros(FP, 1), _ones(FP, 2))
+    buffer = HCubature.hcubature_buffer(f, _zeros(FP, 1), _ones(FP, 1))
     rule = HAdaptiveCubature(rule.kwargs..., buffer = buffer)
 
     # Convert the Rope into Segments, sum the integrals of those

From 1727b79d47f7293a2453f7e1ed09076f77a5060b Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 13:07:10 -0500
Subject: [PATCH 04/15] Use a sample integrand to inform buffer

---
 src/specializations/Rope.jl | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index 594f1a6a..37a2f0ce 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -43,11 +43,17 @@ function integral(
     FP::Type{T} = Float64,
     kwargs...
 ) where {T <: AbstractFloat}
-    # Append a buffer to the given rule
-    buffer = HCubature.hcubature_buffer(f, _zeros(FP, 1), _ones(FP, 1))
+    # Geometry information
+    N = Meshes.paramdim(geometry)
+    segments = Meshes.segments(rope)
+
+    # Use a sample integrand to develop and append a buffer to the given rule
+    integrand(ts) = f(segments[1](ts...)) * differential(segments[1], ts)
+    uintegrand(ts) = Unitful.ustrip.(integrand(ts))
+    buffer = HCubature.hcubature_buffer(uintegrand, _zeros(FP, N), _ones(FP, N))
     rule = HAdaptiveCubature(rule.kwargs..., buffer = buffer)
 
     # Convert the Rope into Segments, sum the integrals of those
     _subintegral(seg) = _integral(f, seg, rule; FP = FP, kwargs...)
-    return sum(_subintegral, Meshes.segments(rope))
+    return sum(_subintegral, segments)
 end

From abdb0ceab9f9753ce52adcf7619a303026d85108 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 13:08:07 -0500
Subject: [PATCH 05/15] Update name

---
 src/specializations/Rope.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index 37a2f0ce..1fb8a8b6 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -36,6 +36,7 @@ function integral(
     return sum(segment -> integral(f, segment, rule; kwargs...), Meshes.segments(rope))
 end
 
+# Use HCubature.hcubature_buffer to reduce allocations
 function integral(
     f,
     rope::Meshes.Rope,
@@ -44,7 +45,7 @@ function integral(
     kwargs...
 ) where {T <: AbstractFloat}
     # Geometry information
-    N = Meshes.paramdim(geometry)
+    N = Meshes.paramdim(rope)
     segments = Meshes.segments(rope)
 
     # Use a sample integrand to develop and append a buffer to the given rule

From bbb6fd6d7eacec6b7193279ff7fc828f9d7b82f4 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 13:45:12 -0500
Subject: [PATCH 06/15] Update accessor to generator

---
 src/specializations/Rope.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index 1fb8a8b6..abd63936 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -49,7 +49,8 @@ function integral(
     segments = Meshes.segments(rope)
 
     # Use a sample integrand to develop and append a buffer to the given rule
-    integrand(ts) = f(segments[1](ts...)) * differential(segments[1], ts)
+    sample = first(segments)
+    integrand(ts) = f(sample(ts...)) * differential(sample, ts)
     uintegrand(ts) = Unitful.ustrip.(integrand(ts))
     buffer = HCubature.hcubature_buffer(uintegrand, _zeros(FP, N), _ones(FP, N))
     rule = HAdaptiveCubature(rule.kwargs..., buffer = buffer)

From b37cff0d8bf64466ee56ccbb0bf7147bdf300226 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 14:16:43 -0500
Subject: [PATCH 07/15] Also benchmark HCubature rules

---
 benchmark/benchmarks.jl | 33 ++++++++++++++++++---------------
 1 file changed, 18 insertions(+), 15 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 14896a03..902d186d 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -41,11 +41,8 @@ end
 spec = (
     f = p -> norm(to(p)),
     f_exp = p::Point -> exp(-norm(to(p))^2 / u"m^2"),
-    g = (
+    geometries = (
         bezier = BezierCurve([Point(t, sin(t), 0) for t in range(-pi, pi, length = 361)]),
-        line = Line(Point(0, 0, 0), Point(1, 1, 1)),
-        plane = Plane(Point(0, 0, 0), Vec(0, 0, 1)),
-        ray = Ray(Point(0, 0, 0), Vec(0, 0, 1)),
         rope = Rope([Point(t, t, t) for t in 1:32]...),
         triangle = Triangle(Point(1, 0, 0), Point(0, 1, 0), Point(0, 0, 1)),
         tetrahedron = let
@@ -56,19 +53,25 @@ spec = (
             Tetrahedron(a, b, c, a + ẑ)
         end
     ),
-    rule_gl = GaussLegendre(100),
-    rule_gk = GaussKronrod(),
-    rule_hc = HAdaptiveCubature()
+    geometries_exp = (
+        line = Line(Point(0, 0, 0), Point(1, 1, 1)),
+        plane = Plane(Point(0, 0, 0), Vec(0, 0, 1)),
+        ray = Ray(Point(0, 0, 0), Vec(0, 0, 1))
+    ),
+    rules = (
+            ( name = "GaussLegendre", rule = GaussLegendre(100) ),
+            ( name = "HAdaptiveCubature", rule = HAdaptiveCubature() )
+    )
 )
 
-SUITE["Specializations/Scalar GaussLegendre"] = let s = BenchmarkGroup()
-    s["BezierCurve"] = @benchmarkable integral($spec.f, $spec.g.bezier, $spec.rule_gl)
-    s["Line"] = @benchmarkable integral($spec.f_exp, $spec.g.line, $spec.rule_gl)
-    s["Plane"] = @benchmarkable integral($spec.f_exp, $spec.g.plane, $spec.rule_gl)
-    s["Ray"] = @benchmarkable integral($spec.f_exp, $spec.g.ray, $spec.rule_gl)
-    s["Rope"] = @benchmarkable integral($spec.f, $spec.g.rope, $spec.rule_gl)
-    s["Triangle"] = @benchmarkable integral($spec.f, $spec.g.triangle, $spec.rule_gl)
-    s["Tetrahedron"] = @benchmarkable integral($spec.f, $spec.g.tetrahedron, $spec.rule_gl)
+SUITE["Specializations"] = let s = BenchmarkGroup()
+    s[]
+    for r in spec.rules, geometry in spec.geometries
+        s[geometry.name, r.name] = @benchmarkable integral($spec.f, geometry, r.rule)
+    end
+    for r in spec.rules, geometry in spec.geometries_exp
+        s[geometry.name, r.name] = @benchmarkable integral($spec.f_exp, geometry, r.rule)
+    end
     s
 end
 

From a04b3ba2c7090b140b897da5e8c09410f020b6bb Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 14:28:00 -0500
Subject: [PATCH 08/15] Generate geometry names

---
 benchmark/benchmarks.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 902d186d..847895f7 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -67,10 +67,12 @@ spec = (
 SUITE["Specializations"] = let s = BenchmarkGroup()
     s[]
     for r in spec.rules, geometry in spec.geometries
-        s[geometry.name, r.name] = @benchmarkable integral($spec.f, geometry, r.rule)
+        geometry_name = nameof(typeof(geometry.name))
+        s[geometry_name, r.name] = @benchmarkable integral($spec.f, geometry, r.rule)
     end
     for r in spec.rules, geometry in spec.geometries_exp
-        s[geometry.name, r.name] = @benchmarkable integral($spec.f_exp, geometry, r.rule)
+        geometry_name = nameof(typeof(geometry.name))
+        s[geometry_name, r.name] = @benchmarkable integral($spec.f_exp, geometry, r.rule)
     end
     s
 end

From b09776ca12527939ac0059ee3195859854d7c738 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 14:34:15 -0500
Subject: [PATCH 09/15] Drop field name

---
 benchmark/benchmarks.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 847895f7..b1d0de52 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -67,11 +67,11 @@ spec = (
 SUITE["Specializations"] = let s = BenchmarkGroup()
     s[]
     for r in spec.rules, geometry in spec.geometries
-        geometry_name = nameof(typeof(geometry.name))
+        geometry_name = nameof(typeof(geometry))
         s[geometry_name, r.name] = @benchmarkable integral($spec.f, geometry, r.rule)
     end
     for r in spec.rules, geometry in spec.geometries_exp
-        geometry_name = nameof(typeof(geometry.name))
+        geometry_name = nameof(typeof(geometry))
         s[geometry_name, r.name] = @benchmarkable integral($spec.f_exp, geometry, r.rule)
     end
     s

From 795be4c37c1a7924b941278f945df0c81d2cb841 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 16:44:17 -0500
Subject: [PATCH 10/15] Try without explicit rule names

---
 benchmark/benchmarks.jl | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index b1d0de52..7f9dfec1 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -59,19 +59,22 @@ spec = (
         ray = Ray(Point(0, 0, 0), Vec(0, 0, 1))
     ),
     rules = (
-            ( name = "GaussLegendre", rule = GaussLegendre(100) ),
-            ( name = "HAdaptiveCubature", rule = HAdaptiveCubature() )
+        GaussLegendre(100),
+        HAdaptiveCubature()
     )
 )
 
 SUITE["Specializations"] = let s = BenchmarkGroup()
-    s[]
+    #=
     for r in spec.rules, geometry in spec.geometries
         geometry_name = nameof(typeof(geometry))
         s[geometry_name, r.name] = @benchmarkable integral($spec.f, geometry, r.rule)
     end
-    for r in spec.rules, geometry in spec.geometries_exp
+    =#
+    for rule in spec.rules, geometry in spec.geometries_exp
+        @info (rule, geometry)
         geometry_name = nameof(typeof(geometry))
+        rule_name = nameof(typeof(rule))
         s[geometry_name, r.name] = @benchmarkable integral($spec.f_exp, geometry, r.rule)
     end
     s

From 2cd7b0c4122e6c6a7919a3e4a241c41b827a3338 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 16:46:58 -0500
Subject: [PATCH 11/15] Update rule name

---
 benchmark/benchmarks.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 7f9dfec1..893d5db9 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -72,10 +72,9 @@ SUITE["Specializations"] = let s = BenchmarkGroup()
     end
     =#
     for rule in spec.rules, geometry in spec.geometries_exp
-        @info (rule, geometry)
         geometry_name = nameof(typeof(geometry))
         rule_name = nameof(typeof(rule))
-        s[geometry_name, r.name] = @benchmarkable integral($spec.f_exp, geometry, r.rule)
+        s[geometry_name, rule_name] = @benchmarkable integral($spec.f_exp, geometry, rule)
     end
     s
 end

From 0f9f1107c24d836bb0ea902916f52b5ff0f74b72 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 18:04:56 -0500
Subject: [PATCH 12/15] Use interp

---
 benchmark/benchmarks.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 893d5db9..dc52337c 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -74,7 +74,7 @@ SUITE["Specializations"] = let s = BenchmarkGroup()
     for rule in spec.rules, geometry in spec.geometries_exp
         geometry_name = nameof(typeof(geometry))
         rule_name = nameof(typeof(rule))
-        s[geometry_name, rule_name] = @benchmarkable integral($spec.f_exp, geometry, rule)
+        s[geometry_name, rule_name] = @benchmarkable integral($spec.f_exp, $geometry, $rule)
     end
     s
 end

From ca1c9ab713071190f89e25bb1e35d3054026d1d3 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 18:24:45 -0500
Subject: [PATCH 13/15] Change dict structures

---
 benchmark/benchmarks.jl | 31 +++++++++++++++----------------
 1 file changed, 15 insertions(+), 16 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index dc52337c..ecf0bcb1 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -16,20 +16,19 @@ integrands = (
     (name = "Vector", f = p -> fill(norm(to(p)), 3))
 )
 rules = (
-    (name = "GaussLegendre", rule = GaussLegendre(100)),
-    (name = "GaussKronrod", rule = GaussKronrod()),
-    (name = "HAdaptiveCubature", rule = HAdaptiveCubature())
+    GaussLegendre(100),
+    GaussKronrod(),
+    HAdaptiveCubature()
 )
 geometries = (
-    (name = "Segment", item = Segment(Point(0, 0, 0), Point(1, 1, 1))),
-    (name = "Sphere", item = Sphere(Point(0, 0, 0), 1.0))
+    Segment(Point(0, 0, 0), Point(1, 1, 1)),
+    Sphere(Point(0, 0, 0), 1.0)
 )
 
 SUITE["Integrals"] = let s = BenchmarkGroup()
     for (int, rule, geometry) in Iterators.product(integrands, rules, geometries)
-        n1 = geometry.name
-        n2 = "$(int.name) $(rule.name)"
-        s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule)
+        name = "$(nameof(typeof(geometry))), $(int.name), $(nameof(typeof(rule)))"
+        s[name] = @benchmarkable integral($int.f, $geometry.item, $rule.rule)
     end
     s
 end
@@ -65,16 +64,16 @@ spec = (
 )
 
 SUITE["Specializations"] = let s = BenchmarkGroup()
-    #=
-    for r in spec.rules, geometry in spec.geometries
-        geometry_name = nameof(typeof(geometry))
-        s[geometry_name, r.name] = @benchmarkable integral($spec.f, geometry, r.rule)
+    # Benchmark most specialization geometries
+    for rule in spec.rules, geometry in spec.geometries
+        name = "$(nameof(typeof(geometry))), Scalar, $(nameof(typeof(rule)))"
+        s[name] = @benchmarkable integral($spec.f, $geometry, $rule)
     end
-    =#
+
+    # Geometries that span an infinite domain use exp function instead
     for rule in spec.rules, geometry in spec.geometries_exp
-        geometry_name = nameof(typeof(geometry))
-        rule_name = nameof(typeof(rule))
-        s[geometry_name, rule_name] = @benchmarkable integral($spec.f_exp, $geometry, $rule)
+        name = "$(nameof(typeof(geometry))), Scalar, $(nameof(typeof(rule)))"
+        s[name] = @benchmarkable integral($spec.f_exp, $geometry, $rule)
     end
     s
 end

From 53e86a8f940af9130c3ee68d464b1537d3919fb5 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 18:27:27 -0500
Subject: [PATCH 14/15] Update naming

---
 benchmark/benchmarks.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index ecf0bcb1..9fe923ac 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -28,7 +28,7 @@ geometries = (
 SUITE["Integrals"] = let s = BenchmarkGroup()
     for (int, rule, geometry) in Iterators.product(integrands, rules, geometries)
         name = "$(nameof(typeof(geometry))), $(int.name), $(nameof(typeof(rule)))"
-        s[name] = @benchmarkable integral($int.f, $geometry.item, $rule.rule)
+        s[name] = @benchmarkable integral($int.f, $geometry, $rule)
     end
     s
 end

From 2e1038d07fe2fd8eb2c5ed4e8162e19ca7be8200 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Thu, 19 Dec 2024 19:58:05 -0500
Subject: [PATCH 15/15] Implement for CylinderSurface

---
 benchmark/benchmarks.jl                |  3 +++
 src/specializations/CylinderSurface.jl | 30 ++++++++++++++++++++++++++
 2 files changed, 33 insertions(+)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 9fe923ac..6a750efa 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -21,6 +21,9 @@ rules = (
     HAdaptiveCubature()
 )
 geometries = (
+    let ẑ = Vec(0, 0, 1)
+        CylinderSurface(Plane(Point(0, 0, 0), ẑ), Plane(Point(0, 0, 3), ẑ), 2.5)
+    end,
     Segment(Point(0, 0, 0), Point(1, 1, 1)),
     Sphere(Point(0, 0, 0), 1.0)
 )
diff --git a/src/specializations/CylinderSurface.jl b/src/specializations/CylinderSurface.jl
index c359e3af..f6375edd 100644
--- a/src/specializations/CylinderSurface.jl
+++ b/src/specializations/CylinderSurface.jl
@@ -31,3 +31,33 @@ function integral(
 
     return sides + top + bottom
 end
+
+function integral(
+    f,
+    cyl::Meshes.CylinderSurface,
+    rule::HAdaptiveCubature;
+    FP::Type{T} = Float64,
+    diff_method::DM = _default_diff_method(cyl, FP),
+    kwargs...
+) where {T <: AbstractFloat, DM <: DifferentiationMethod}
+    _check_diff_method_support(cyl, diff_method)
+
+    # Use a sample integrand to develop and append a buffer to the given rule
+    f_sample(ts) = Unitful.ustrip.(f(cyl(ts...)) * differential(cyl, ts))
+    N = Meshes.paramdim(cyl)
+    buffer = HCubature.hcubature_buffer(f_sample, _zeros(FP, N), _ones(FP, N))
+    rule = HAdaptiveCubature(rule.kwargs..., buffer = buffer)
+
+    # The generic method only parametrizes the sides
+    sides = _integral(f, cyl, rule; diff_method = diff_method, FP = FP, kwargs...)
+
+    # Integrate the Disk at the top
+    disk_top = Meshes.Disk(cyl.top, cyl.radius)
+    top = _integral(f, disk_top, rule; diff_method = diff_method, FP = FP, kwargs...)
+
+    # Integrate the Disk at the bottom
+    disk_bottom = Meshes.Disk(cyl.bot, cyl.radius)
+    bottom = _integral(f, disk_bottom, rule; diff_method = diff_method, FP = FP, kwargs...)
+
+    return sides + top + bottom
+end