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] 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