Updated conic integration methods

Project.toml

@@ -1,7 +1,7 @@
 name = "MeshIntegrals"
 uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
 authors = ["Mike Ingold <[email protected]>"]
-version = "0.12.2"
+version = "0.13.0"
 
 [deps]
 CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
@@ -13,11 +13,11 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
-CoordRefSystems = "0.10, 0.11"
+CoordRefSystems = "0.12, 0.13"
 FastGaussQuadrature = "1"
 HCubature = "1.5"
 LinearAlgebra = "1"
-Meshes = "0.47 - 0.49"
+Meshes = "0.50"
 QuadGK = "2"
 Unitful = "1"
 julia = "1.6"

README.md

@@ -1,10 +1,13 @@
 # MeshIntegrals.jl
 
+[![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT)
+[![ColPrac](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet?style=flat-square)](https://github.com/SciML/ColPrac)
+
 [![Build Status](https://github.com/mikeingold/MeshIntegrals.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/mikeingold/MeshIntegrals.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![codecov](https://codecov.io/gh/mikeingold/MeshIntegrals.jl/graph/badge.svg)](https://codecov.io/gh/mikeingold/MeshIntegrals.jl)
 [![Coveralls](https://coveralls.io/repos/github/mikeingold/MeshIntegrals.jl/badge.svg?branch=main)](https://coveralls.io/github/mikeingold/MeshIntegrals.jl?branch=main)
 [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
-[![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT)
+
 
 This package implements methods for numerically-computing integrals over geometric polytopes
 from [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) using the following `::IntegrationAlgorithms`:
@@ -51,8 +54,7 @@ integral(f, unit_circle_bz, GaussKronrod())
 | Symbol | Meaning |
 |--------|---------|
 | :white_check_mark: | Implemented, passes tests |
-| :x: | Planned but not yet implemented |
-| :warning: | Unable to implement: parameterization not available (see [Issue #28](https://github.com/mikeingold/MeshIntegrals.jl/issues/28)) |
+| :x: | Not yet supported |
 | :stop_sign: | Not supported |
 
 ### Integral
@@ -64,15 +66,14 @@ integral(f, unit_circle_bz, GaussKronrod())
 | `Box` in `𝔼{1}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Box` in `𝔼{2}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Box` in `𝔼{3}` | :white_check_mark: | :stop_sign: | :white_check_mark: |
-| `Box` in `𝔼{N}` | :x: | :x: | :x: |
 | `Circle` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
-| `Cone` | :warning: | :warning: | :warning: |
-| `ConeSurface` | :x: | :x: | :x: |
+| `Cone` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
+| `ConeSurface` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Cylinder` | :white_check_mark: | :stop_sign: | :white_check_mark: |
-| `CylinderSurface` | :x: | :white_check_mark: | :white_check_mark: |
+| `CylinderSurface` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Disk` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
-| `Frustum` | :warning: | :warning: | :warning: |
-| `FrustumSurface` | :x: | :x: | :x: |
+| `Frustum` | :stop_sign: | :stop_sign: | :stop_sign: |
+| `FrustumSurface` | :stop_sign: | :stop_sign: | :stop_sign: |
 | `Line` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `ParaboloidSurface` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Plane` | :white_check_mark: | :white_check_mark: | :white_check_mark: |

src/integral_surface.jl

@@ -33,91 +33,109 @@
     return QuadGK.quadgk(v -> ∫₁(v), FP(0), FP(1); settings.kwargs...)[1]
 end
 
-function _integral_2d(
-    f,
-    geometry,
-    settings::HAdaptiveCubature,
-)
-    return _integral(f, geometry, settings)
-end
-
-function _integral_2d(
-    f,
-    geometry,
-    settings::HAdaptiveCubature,
-    FP::Type{T}
-) where {T<:AbstractFloat}
-    return _integral(f, geometry, settings, FP)
-end
-
 
 ################################################################################
-#                  Specialized Methods for CylinderSurface
+#    Specialized Methods for ConeSurface, CylinderSurface, and FrustumSurface
 ################################################################################
 
 function integral(
     f::F,
     cyl::Meshes.CylinderSurface,
-    settings::GaussLegendre,
+    settings::I,
     FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
-    error("Integrating a CylinderSurface with GaussLegendre not supported.")
-    # TODO Planned to support in the future
+) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
+    # The generic method only parameterizes the sides
+    sides = _integral_2d(f, cyl, settings, FP)
+
+    # Integrate the Disk at the top
+    disk_top = Meshes.Disk(cyl.top, cyl.radius)
+    top = _integral_2d(f, disk_top, settings, FP)
+
+    # Integrate the Disk at the bottom
+    disk_bottom = Meshes.Disk(cyl.bot, cyl.radius)
+    bottom = _integral_2d(f, disk_bottom, settings, FP)
+
+    return sides + top + bottom
 end
 
 function integral(
     f::F,
     cyl::Meshes.CylinderSurface,
-    settings::GaussKronrod,
+    settings::HAdaptiveCubature,
     FP::Type{T} = Float64
 ) where {F<:Function, T<:AbstractFloat}
-    # Integrate the curved sides of the CylinderSurface
-    # \int ( \int f(r̄) dz ) dφ
-    function sides_inner∫(φ)
-        sidelength = norm(cyl(φ,FP(1)) - cyl(φ,FP(0)))
-        return sidelength * QuadGK.quadgk(z -> f(cyl(φ,z)), FP(0), FP(1); settings.kwargs...)[1]
-    end
-    sides = (FP(2π) * cyl.radius) .* QuadGK.quadgk(φ -> sides_inner∫(φ), FP(0), FP(1); settings.kwargs...)[1]
+    # The generic method only parameterizes the sides
+    sides = _integral(f, cyl, settings, FP)
 
-    # Integrate the Disk at the top of the CylinderSurface
+    # Integrate the Disk at the top
     disk_top = Meshes.Disk(cyl.top, cyl.radius)
-    top = _integral_2d(f, disk_top, settings, FP)
+    top = _integral(f, disk_top, settings, FP)
 
-    # Integrate the Disk at the bottom of the CylinderSurface
+    # Integrate the Disk at the bottom
     disk_bottom = Meshes.Disk(cyl.bot, cyl.radius)
-    bottom = _integral_2d(f, disk_bottom, settings, FP)
+    bottom = _integral(f, disk_bottom, settings, FP)
 
     return sides + top + bottom
 end
 
 function integral(
     f::F,
-    cyl::Meshes.CylinderSurface,
+    cone::Meshes.ConeSurface,
+    settings::I,
+    FP::Type{T} = Float64
+) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
+    # The generic method only parameterizes the sides
+    sides = _integral_2d(f, cone, settings, FP)
+
+    # Integrate the Disk at the base
+    base = _integral_2d(f, cone.base, settings, FP)
+
+    return sides + base
+end
+
+function integral(
+    f::F,
+    cone::Meshes.ConeSurface,
     settings::HAdaptiveCubature,
     FP::Type{T} = Float64
 ) where {F<:Function, T<:AbstractFloat}
-    Dim = Meshes.paramdim(cyl)
+    # The generic method only parameterizes the sides
+    sides = _integral(f, cone, settings, FP)
 
-    integrand(t) = f(cyl(t...)) * differential(cyl, t)
+    # Integrate the Disk at the base
+    base = _integral(f, cone.base, settings, FP)
 
-    # HCubature doesn't support functions that output Unitful Quantity types
-    # Establish the units that are output by f
-    testpoint_parametriccoord = fill(FP(0.5),Dim)
-    integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
-    # Create a wrapper that returns only the value component in those units
-    uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
-    # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, zeros(FP,Dim), ones(FP,Dim); settings.kwargs...)[1]
-    # Reapply units
-    sides = value .* integrandunits
+    return sides + base
+end
 
-    # Integrate the Disk at the top of the CylinderSurface
-    disk_top = Meshes.Disk(cyl.top, cyl.radius)
-    top = _integral_2d(f, disk_top, settings, FP)
+function integral(
+    f::F,
+    frust::Meshes.FrustumSurface,
+    settings::I,
+    FP::Type{T} = Float64
+) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
+    # The generic method only parameterizes the sides
+    sides = _integral_2d(f, frust, settings, FP)
 
-    # Integrate the Disk at the bottom of the CylinderSurface
-    disk_bottom = Meshes.Disk(cyl.bot, cyl.radius)
-    bottom = _integral_2d(f, disk_bottom, settings, FP)
+    # Integrate the Disks at the top and bottom
+    top = _integral_2d(f, frust.top, settings, FP)
+    bottom = _integral_2d(f, frust.bot, settings, FP)
+
+    return sides + top + bottom
+end
+
+function integral(
+    f::F,
+    frust::Meshes.FrustumSurface,
+    settings::HAdaptiveCubature,
+    FP::Type{T} = Float64
+) where {F<:Function, T<:AbstractFloat}
+    # The generic method only parameterizes the sides
+    sides = _integral(f, frust, settings, FP)
+
+    # Integrate the Disks at the top and bottom
+    top = _integral(f, frust.top, settings, FP)
+    bottom = _integral(f, frust.bot, settings, FP)
 
     return sides + top + bottom
 end

src/integral_volume.jl

@@ -24,21 +24,14 @@ function _integral_3d(
     return FP(1/8) .* sum(integrand, zip(wws,xxs))
 end
 
+# Integrating volumes with GaussKronrod not supported by default
 function _integral_3d(
     f,
     geometry,
-    settings::HAdaptiveCubature,
-)
-    return _integral(f, geometry, settings)
-end
-
-function _integral_3d(
-    f,
-    geometry,
-    settings::HAdaptiveCubature,
-    FP::Type{T}
+    settings::GaussKronrod,
+    FP::Type{T} = Float64
 ) where {T<:AbstractFloat}
-    return _integral(f, geometry, settings, FP)
+    error("Integrating this volume type with GaussKronrod not supported.")
 end
 
 
@@ -77,35 +70,3 @@ function integral(
 ) where {F<:Function, T<:AbstractFloat}
     error("Integrating a Tetrahedron with HAdaptiveCubature not supported.")
 end
-
-
-################################################################################
-#                         Unsupported Placeholders
-################################################################################
-
-function integral(
-    f::F,
-    ball::Meshes.Ball{Meshes.𝔼{3},CRS,ℒ},
-    settings::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, CRS, ℒ, T<:AbstractFloat}
-    error("Integrating a Ball in 𝔼{3} with GaussKronrod not supported.")
-end
-
-function integral(
-    f::F,
-    box::Meshes.Box{Meshes.𝔼{3},CRS},
-    settings::GaussKronrod,
-    FP::Type{T} = Float64,
-) where {F<:Function, CRS, T<:AbstractFloat}
-    error("Integrating a Box in 𝔼{3} with GaussKronrod not supported.")
-end
-
-function integral(
-    f::F,
-    box::Meshes.Cylinder,
-    settings::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
-    error("Integrating a Cylinder with GaussKronrod not supported.")
-end

test/Project.toml

@@ -8,7 +8,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 Aqua = "0.8"
-Meshes = "0.49"
-MeshIntegrals = "0.12"
+Meshes = "0.50"
+MeshIntegrals = "0.13"
 QuadGK = "2"
 Unitful = "1"