Add more integral methods

* Add surfaceintegral method for CylinderSurface using GaussKronrod

* Add tests for surfaceintegral on CylinderSurface

* Update support matrix

* Add temporary debugging

* Bugfix in math

* Adjust tests

* Remove debugging

* Add volumeintegral method for Ball{3,T}

* Disable complex test due to change in Meshes

* Update support matrix

* Add more surfaceintegral methods for Sphere{3,T}

* Add another volumeintegral method for Ball{3,T}

* Bugfix

* Move function to another section

* Bump version and update support matrix

* Remove branch-specific plans

Author: mikeingold

Project.toml

@@ -1,7 +1,7 @@
 name = "MeshIntegrals"
 uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
 authors = ["Mike Ingold <[email protected]>"]
-version = "0.9.4"
+version = "0.9.5"
 
 [deps]
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"

README.md

@@ -20,10 +20,15 @@ Methods are tested to ensure compatibility with
 | Symbol | Meaning |
 |--------|---------|
 | :white_check_mark: | Implemented, passes tests |
-| :yellow_square: | Implemented, untested |
+| :yellow_square: | Implemented, not yet validated |
 | :x: | Not yet implemented |
 
-### Line Integrals
+### Integral
+| Geometry | Gauss-Legendre | Gauss-Kronrod | H-Adaptive Cubature |
+|----------|----------------|---------------|---------------------|
+| `Meshes.Box{Dim,T}` | :x: | :x: | :x: |
+
+### Line Integral
 | Geometry | Gauss-Legendre | Gauss-Kronrod | H-Adaptive Cubature |
 |----------|----------------|---------------|---------------------|
 | `Meshes.BezierCurve` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
@@ -35,25 +40,24 @@ Methods are tested to ensure compatibility with
 | `Meshes.Segment` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Meshes.Sphere{2,T}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 
-### Surface Integrals
+### Surface Integral
 | Geometry | Gauss-Legendre | Gauss-Kronrod | H-Adaptive Cubature |
 |----------|----------------|---------------|-------------------|
 | `Meshes.Ball{2,T}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Meshes.Box{2,T}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
-| `Meshes.CylinderSurface` | :x: | :x: | :x: |
+| `Meshes.CylinderSurface` | :x: | :white_check_mark: | :x: |
 | `Meshes.Disk` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Meshes.ParaboloidSurface` | :x: | :x: | :x: |
-| `Meshes.Sphere{3,T}` | :x: | :x: | :x: |
+| `Meshes.Sphere{3,T}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Meshes.Triangle` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Meshes.Torus` | :x: | :x: | :x: |
 | `Meshes.SimpleMesh` | :x: | :x: | :x: |
 
-### Volume Integrals
+### Volume Integral
 | Geometry | Gauss-Legendre | H-Adaptive Cubature |
 |----------|----------------|---------------|
-| `Meshes.Ball{3,T}` | :x: | :x: |
+| `Meshes.Ball{3,T}` | :white_check_mark: | :white_check_mark: |
 | `Meshes.Box{3,T}` | :white_check_mark: | :white_check_mark: |
-| `Meshes.Box{Dim,T}` | :x: | :x: |
 
 # Example Usage
 

src/integral_surface.jl

@@ -127,6 +127,30 @@ function surfaceintegral(
     return (π/4) * area(triangle) .* sum(integrand, zip(wws,xxs))
 end
 
+function surfaceintegral(
+    f,
+    sphere::Meshes.Sphere{3,T},
+    settings::GaussLegendre
+) where {T}
+    # Validate the provided integrand function
+    _validate_integrand(f,3,T)
+
+    # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2
+    xs, ws = gausslegendre(settings.n)
+    wws = Iterators.product(ws, ws)
+    xxs = Iterators.product(xs, xs)
+
+    # Domain transformation: xi,xj [-1,1] ↦ s,t [0,1]
+    t(xi) = 0.5xi + 0.5
+    u(xj) = 0.5xj + 0.5
+
+    # Integrate the sphere in parametric (t,u)-space [0,1]²
+    integrand(t,u) = sinpi(t) * f(sphere(1,t,u))
+    g(((wi,wj), (xi,xj))) = wi * wj * integrand(t(xi),u(xj))
+    R = sphere.radius
+    return 0.25 * 2π^2 * R^2 .* sum(g, zip(wws,xxs))
+end
+
 ################################################################################
 #                               Gauss-Kronrod
 ################################################################################
@@ -167,6 +191,41 @@ function surfaceintegral(
     return area(box) .* outerintegral
 end
 
+function surfaceintegral(
+    f,
+    cyl::Meshes.CylinderSurface{T},
+    settings::GaussKronrod
+) where {T}
+    # Validate the provided integrand function
+    # A CylinderSurface is definitionally embedded in 3D-space
+    _validate_integrand(f,3,T)
+
+    # Integrate the rounded sides of the cylinder's surface
+    # \int ( \int f(r̄) dz ) dφ
+    function sides_innerintegral(φ)
+        sidelength = norm(cyl(φ,1) - cyl(φ,0))
+        return sidelength * quadgk(z -> f(cyl(φ,z)), 0, 1; settings.kwargs...)[1]
+    end
+    sides = (2π * cyl.radius) .* quadgk(φ -> sides_innerintegral(φ), 0, 1; settings.kwargs...)[1]
+
+    # Integrate the top and bottom disks
+    # \int ( \int r f(r̄) dr ) dφ
+    function disk_innerintegral(φ,plane,z)
+        # Parameterize the top surface of the cylinder
+        rimedge = cyl(φ,z)
+        centerpoint = plane.p
+        r̄ = rimedge - centerpoint
+        radius = norm(r̄)
+        point(r) = centerpoint + (r / radius) * r̄
+
+        return radius^2 * quadgk(r -> r * f(point(r)), 0, 1; settings.kwargs...)[1]
+    end
+    top    = 2π .* quadgk(φ -> disk_innerintegral(φ,cyl.top,1), 0, 1; settings.kwargs...)[1]
+    bottom = 2π .* quadgk(φ -> disk_innerintegral(φ,cyl.bot,0), 0, 1; settings.kwargs...)[1]
+
+    return sides + top + bottom
+end
+
 function surfaceintegral(
     f,
     disk::Meshes.Disk{T},
@@ -212,6 +271,22 @@ function surfaceintegral(
     return 2.0 * area(triangle) .* outerintegral
 end
 
+function surfaceintegral(
+    f,
+    sphere::Meshes.Sphere{3,T},
+    settings::GaussKronrod
+) where {T}
+    # Validate the provided integrand function
+    _validate_integrand(f,3,T)
+
+    # Integrate the sphere in parametric (t,u)-space [0,1]^2
+    innerintegrand(u) = quadgk(t -> sinpi(t) * f(sphere(1,t,u)), 0, 1)[1]
+    intval = quadgk(u -> innerintegrand(u), 0, 1, settings.kwargs...)[1]
+
+    R = sphere.radius
+    return 2π^2 * R^2 .* intval
+end
+
 
 ################################################################################
 #                               HCubature
@@ -307,3 +382,20 @@ function surfaceintegral(
     # Apply a linear domain-correction factor 0.5 ↦ area(triangle)
     return 2.0 * area(triangle) .* intval
 end
+
+function surfaceintegral(
+    f,
+    sphere::Meshes.Sphere{3,T},
+    settings::HAdaptiveCubature
+) where {T}
+    # Validate the provided integrand function
+    _validate_integrand(f,3,T)
+
+    # Integrate the sphere in parametric (t,u)-space [0,1]^2
+    integrand(t,u) = sinpi(t) * f(sphere(1,t,u))
+    integrand(tu) = integrand(tu[1],tu[2])
+    intval = hcubature(tu -> integrand(tu), [0,0], [1,1], settings.kwargs...)[1]
+
+    R = sphere.radius
+    return 2π^2 * R^2 .* intval
+end

src/integral_volume.jl

@@ -2,6 +2,31 @@
 #                               Gauss-Legendre
 ################################################################################
 
+function volumeintegral(
+    f,
+    ball::Meshes.Ball{3,T},
+    settings::GaussLegendre
+) where {T}
+    # Validate the provided integrand function
+    _validate_integrand(f,3,T)
+
+    # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2
+    xs, ws = gausslegendre(settings.n)
+    wws = Iterators.product(ws, ws, ws)
+    xxs = Iterators.product(xs, xs, xs)
+
+    # Domain transformation: xi,xj,xk [-1,1] ↦ s,t,u [0,1]
+    s(xi) = 0.5xi + 0.5
+    t(xj) = 0.5xj + 0.5
+    u(xk) = 0.5xk + 0.5
+
+    # Integrate the ball in parametric (s,t,u)-space [0,1]³
+    integrand(s,t,u) = s^2 * sinpi(t) * f(ball(s,t,u))
+    g(((wi,wj,wk), (xi,xj,xk))) = wi * wj * wk * integrand(s(xi),t(xj),u(xk))
+    R = ball.radius
+    return (1/8) * 2π^2 * R^3 .* sum(g, zip(wws,xxs))
+end
+
 function volumeintegral(
     f,
     box::Meshes.Box{3,T},
@@ -34,6 +59,23 @@ end
 #                               HCubature
 ################################################################################
 
+function volumeintegral(
+    f,
+    ball::Meshes.Ball{3,T},
+    settings::HAdaptiveCubature
+) where {T}
+    # Validate the provided integrand function
+    _validate_integrand(f,3,T)
+
+    # Integrate the ball in parametric (s,t,u)-space [0,1]³
+    integrand(s,t,u) = s^2 * sinpi(t) * f(ball(s,t,u))
+    integrand(stu) = integrand(stu[1],stu[2],stu[3])
+    intval = hcubature(stu -> integrand(stu), [0,0,0], [1,1,1], settings.kwargs...)[1]
+
+    R = ball.radius
+    return 2π^2 * R^3 .* intval
+end
+
 function volumeintegral(
     f,
     box::Meshes.Box{3,T},

test/runtests.jl

@@ -23,9 +23,6 @@ using Test
 
     # Line segments/paths oriented CCW between points
     seg_ne = Segment(pt_e, pt_n)
-    #seg_nw = Segment(pt_n, pt_w)
-    #seg_sw = Segment(pt_w, pt_s)
-    #seg_se = Segment(pt_s, pt_e)
     line_ne = Line(pt_e, pt_n)
     ring_rect = Ring(pt_e, pt_n, pt_w, pt_s)
     rope_rect = Rope(pt_e, pt_n, pt_w, pt_s, pt_e)
@@ -44,6 +41,8 @@ using Test
     sphere2d = Sphere(origin2d, 2.0)
     sphere3d = Sphere(origin3d, 2.0)
     triangle = Ngon(pt_e, pt_n, pt_w)
+    cylsurf = CylinderSurface(pt_e, pt_w, 2.5)
+    # TODO add test for a non-right-cylinder surface when measure(c) is fixed in Meshes
 
     @testset "Errors Expected on Invalid Methods" begin
         f(::Point) = 1.0
@@ -53,7 +52,7 @@ using Test
         @test_throws MethodError lineintegral(f, ball3d)      # Ball{3,T}
         @test_throws MethodError lineintegral(f, box2d)       # Box{2,T}
         @test_throws MethodError lineintegral(f, box3d)       # Box{3,T}
-        # CylinderSurface
+        @test_throws MethodError lineintegral(f, cylsurf)     # CylinderSurface
         @test_throws MethodError lineintegral(f, disk)        # Disk
         @test_throws MethodError lineintegral(f, triangle)    # Ngon{3,Dim,T}
         # ParaboloidSurface
@@ -101,9 +100,11 @@ using Test
                 @test isapprox(surfaceintegral(f, box2d, rule), area(box2d); rtol=1e-6)         # Box{2,T}
                 @test isapprox(surfaceintegral(f, disk, rule), area(disk); rtol=1e-6)           # Disk
                 @test isapprox(surfaceintegral(f, triangle, rule), area(triangle); rtol=1e-6)   # Triangle
+                @test isapprox(surfaceintegral(f, cylsurf, rule), area(cylsurf); rtol=1e-6)     # CylinderSurface
 
                 # Volume Integrals (skip for GaussKronrod rules)
                 if rule != GaussKronrod()
+                    @test volumeintegral(f, ball3d, rule) ≈ volume(ball3d)   # Ball{3,T}
                     @test volumeintegral(f, box3d, rule) ≈ volume(box3d)     # Box{3,T}
                 end
             end
@@ -127,15 +128,18 @@ using Test
                 @test isapprox(surfaceintegral(f, box2d, rule), fill(area(box2d),3); rtol=1e-6)        # Box{2,T}
                 @test isapprox(surfaceintegral(f, disk, rule), fill(area(disk),3); rtol=1e-6)          # Disk
                 @test isapprox(surfaceintegral(f, triangle, rule), fill(area(triangle),3); rtol=1e-6)  # Triangle
+                @test isapprox(surfaceintegral(f, cylsurf, rule), fill(area(cylsurf),3); rtol=1e-6)    # CylinderSurface
 
                 # Volume Integrals (skip for GaussKronrod rules)
                 if rule != GaussKronrod()
+                    @test volumeintegral(f, ball3d, rule) ≈ fill(volume(ball3d),3)   # Ball{3,T}
                     @test volumeintegral(f, box3d, rule) ≈ fill(volume(box3d),3)     # Box{3,T}
                 end
             end
         end
     end
 
+    #= Disabled: As of Feb 2024 Meshes seems to now disallow Point{1,Complex}
     @testset "Contour Integrals on a Point{1,Complex}-Domain" begin
         fc(z::Complex) = 1/z
         fc(p::Point{1,ComplexF64}) = fc(p.coords[1])
@@ -148,4 +152,5 @@ using Test
         # ∴ \int_C (1/z) dz = 2πi
         @test lineintegral(fc, unit_circle_complex, GaussKronrod()) ≈ 2π*im
     end
+    =#
 end