Obsolete nested GaussKronrod rules (#196)

* Remove support for GK in more than one dimension

* Bump version to 17-dev

* Update docs and changelog

* Update support status for surfaces

* Remove unused deprecated status entry from legend

* Tweak error check

* Tweak other error checks

* Update benchmarks to skip unsupported GK usages

* Update outdated entries

* Use helper function for alias error message generation

* Revert error message

Author: mikeingold

CHANGELOG.md

@@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Removed
+
+- Removed previously-deprecated support for use of `GaussKronrod` rules on geometries with more than one parametric dimension.
+
 
 ## [0.16.4] - 2025-08-10
 

Project.toml

@@ -1,6 +1,6 @@
 name = "MeshIntegrals"
 uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
-version = "0.16.4"
+version = "0.17.0-DEV"
 
 [deps]
 CliffordNumbers = "3998ac73-6bd4-4031-8035-f167dd3ed523"

benchmark/benchmarks.jl

@@ -21,12 +21,19 @@ rules = (
     (name = "HAdaptiveCubature", rule = HAdaptiveCubature())
 )
 geometries = (
-    (name = "Segment", item = Segment(Point(0, 0, 0), Point(1, 1, 1))),
-    (name = "Sphere", item = Sphere(Point(0, 0, 0), 1.0))
+    (name = "Segment", item = Segment(Point(0, 0, 0), Point(1, 1, 1))), # 1D
+    (name = "Sphere", item = Sphere(Point(0, 0), 1.0)), # 2D
+    (name = "Sphere", item = Sphere(Point(0, 0, 0), 1.0)) # 3D
 )
 
 SUITE["Integrals"] = let s = BenchmarkGroup()
     for (int, rule, geometry) in Iterators.product(integrands, rules, geometries)
+        # Skip unsupported applications of GaussKronrod
+        if (Meshes.paramdim(geometry.item) > 1) && (rule.rule == GaussKronrod())
+            continue
+        end
+
+        # Construct benchmark and add it to test suite
         n1 = geometry.name
         n2 = "$(int.name) $(rule.name)"
         s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule)

docs/src/support.md

@@ -2,24 +2,13 @@
 
 This library aims to enable users to calculate the value of integrals over all
 [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) geometry types using
-a number of numerical integration rules and techniques. However, some combinations
-of geometry types and integration rules are ill-suited (and a few are simply not
-yet implemented).
-
-## General Recommendations
+a number of numerical integration rules and techniques.
 
 In general, `GaussKronrod` integration rules are recommended (and the default) for
 geometries with one parametric dimension. For geometries with more than one
 parametric dimension, e.g. surfaces and volumes, `HAdaptiveCubature` rules are
 recommended (and the default).
 
-While it is currently possible to apply nested `GaussKronrod` rules to numerically
-integrate surfaces, this produces results that are strictly inferior to using an
-equivalent `HAdaptiveCubature` rule, so support for this usage has been deprecated.
-In version 16.x of MeshIntegrals.jl, using a `GaussKronrod` rule for a surface
-will work but will yield a deprecation warning. Beginning with a future version
-17.0, this combination will simply be unsupported and throw an error.
-
 ## The Support Matrix
 
 The following Support Matrix captures the current state of support for all geometry/rule
@@ -28,46 +17,47 @@ designed to check for accuracy.
 
 | `Meshes.Geometry/Domain` | `GaussKronrod` | `GaussLegendre` | `HAdaptiveCubature` |
 |----------|----------------|---------------|---------------------|
-| `Ball` in `𝔼{2}` | ⚠️ | ✅ | ✅ |
+| `Ball` in `𝔼{2}` | 🛑 | ✅ | ✅ |
 | `Ball` in `𝔼{3}` | 🛑 | ✅ | ✅ |
 | `BezierCurve` | ✅ | ✅ | ✅ |
 | `Box` in `𝔼{1}` | ✅ | ✅ | ✅ |
-| `Box` in `𝔼{2}` | ⚠️ | ✅ | ✅ |
+| `Box` in `𝔼{2}` | 🛑 | ✅ | ✅ |
 | `Box` in `𝔼{≥3}` | 🛑 | ✅ | ✅ |
-| `CartesianGrid` | ✅ | ✅ | ✅ |
+| `CartesianGrid` in `𝔼{1}` | ✅ | ✅ | ✅ |
+| `CartesianGrid` in `𝔼{≥2}` | 🛑 | ✅ | ✅ |
 | `Circle` | ✅ | ✅ | ✅ |
 | `Cone` | 🛑 | ✅ | ✅ |
-| `ConeSurface` | ⚠️ | ✅ | ✅ |
+| `ConeSurface` | 🛑 | ✅ | ✅ |
 | `Cylinder` | 🛑 | ✅ | ✅ |
-| `CylinderSurface` | ⚠️ | ✅ | ✅ |
-| `Disk` | ⚠️ | ✅ | ✅ |
-| `Ellipsoid` | ✅ | ✅ | ✅ |
-| `Frustum` | ⚠️ | ✅ | ✅ |
-| `FrustumSurface` | ⚠️ | ✅ | ✅ |
-| `Hexahedron` | ✅ | ✅ | ✅ |
+| `CylinderSurface` | 🛑 | ✅ | ✅ |
+| `Disk` | 🛑 | ✅ | ✅ |
+| `Ellipsoid` | 🛑 | ✅ | ✅ |
+| `Frustum` | 🛑 | ✅ | ✅ |
+| `FrustumSurface` | 🛑 | ✅ | ✅ |
+| `Hexahedron` | 🛑 | ✅ | ✅ |
 | `Line` | ✅ | ✅ | ✅ |
-| `ParaboloidSurface` | ⚠️ | ✅ | ✅ |
+| `ParaboloidSurface` | 🛑 | ✅ | ✅ |
 | `ParametrizedCurve` | ✅ | ✅ | ✅ |
-| `Plane` | ✅ | ✅ | ✅ |
-| `PolyArea` | ⚠️ | ✅ | ✅ |
-| `Pyramid` | ⚠️ | ✅ | ✅ |
-| `Quadrangle` | ⚠️ | ✅ | ✅ |
+| `Plane` | 🛑 | ✅ | ✅ |
+| `PolyArea` | 🛑 | ✅ | ✅ |
+| `Pyramid` | 🛑 | ✅ | ✅ |
+| `Quadrangle` | 🛑 | ✅ | ✅ |
 | `Ray` | ✅ | ✅ | ✅ |
-| `RegularGrid` | ✅ | ✅ | ✅ |
+| `RegularGrid` in `𝔼{1}` | ✅ | ✅ | ✅ |
+| `RegularGrid` in `𝔼{≥2}` | 🛑 | ✅ | ✅ |
 | `Ring` | ✅ | ✅ | ✅ |
 | `Rope` | ✅ | ✅ | ✅ |
 | `Segment` | ✅ | ✅ | ✅ |
-| `SimpleMesh` | ⚠️ | ✅ | ✅ |
-| `Sphere` in `𝔼{2}` | ✅ | ✅ | ✅ |
-| `Sphere` in `𝔼{3}` | ⚠️ | ✅ | ✅ |
-| `StructuredGrid` | ✅ | ✅ | ✅ |
-| `Tetrahedron` | ⚠️ | ✅ | ✅ |
-| `Triangle` | ✅ | ✅ | ✅ |
-| `Torus` | ⚠️ | ✅ | ✅ |
-| `Wedge` | ⚠️ | ✅ | ✅ |
+| `SimpleMesh` | 🛑 | ✅ | ✅ |
+| `Sphere` | 🛑 | ✅ | ✅ |
+| `StructuredGrid` in `𝔼{1}` | ✅ | ✅ | ✅ |
+| `StructuredGrid` in `𝔼{≥2}` | 🛑 | ✅ | ✅ |
+| `Tetrahedron` | 🛑 | ✅ | ✅ |
+| `Triangle` | 🛑 | ✅ | ✅ |
+| `Torus` | 🛑 | ✅ | ✅ |
+| `Wedge` | 🛑 | ✅ | ✅ |
 
 | Symbol | Support Level |
 |--------|---------|
 | ✅ | Supported |
-| ⚠️ | Deprecated |
 | 🛑 | Not supported |

src/integral.jl

@@ -63,24 +63,14 @@ function _integral(
 ) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     _check_diff_method_support(geometry, diff_method)
 
-    # Implementation depends on number of parametric dimensions over which to integrate
-    N = Meshes.paramdim(geometry)
-    if N == 1
-        integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method)
-        return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
-    elseif N == 2
-        # Issue deprecation warning
-        Base.depwarn("Use `HAdaptiveCubature` instead of \
-                     `GaussKronrod` for surfaces.", :integral)
-
-        # Nested integration
-        integrand2D(u, v) = f(geometry(u, v)) * differential(geometry, (u, v), diff_method)
-        ∫₁(v) = QuadGK.quadgk(u -> integrand2D(u, v), zero(FP), one(FP); rule.kwargs...)[1]
-        return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1]
-    else
-        _error_unsupported_combination("geometry with more than two parametric dimensions",
-            "GaussKronrod")
+    # Only supported for 1D geometries
+    if Meshes.paramdim(geometry) != 1
+        msg = "GaussKronrod rules not supported in more than one parametric dimension."
+        throw(ArgumentError(msg))
     end
+
+    integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method)
+    return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
 end
 
 # GaussLegendre

src/integral_aliases.jl

@@ -1,3 +1,7 @@
+function _alias_error_msg(name, N)
+    "Performing a $name integral on a geometry with $N parametric dimensions not supported."
+end
+
 ################################################################################
 #                              Line Integral
 ################################################################################
@@ -24,8 +28,7 @@ function lineintegral(
         kwargs...
 )
     if (N = Meshes.paramdim(geometry)) != 1
-        throw(ArgumentError("Performing a line integral on a geometry \
-                            with $N parametric dimensions not supported."))
+        throw(ArgumentError(_alias_error_msg("line", N)))
     end
 
     return integral(f, geometry, rule; kwargs...)
@@ -57,8 +60,7 @@ function surfaceintegral(
         kwargs...
 )
     if (N = Meshes.paramdim(geometry)) != 2
-        throw(ArgumentError("Performing a surface integral on a geometry \
-                            with $N parametric dimensions not supported."))
+        throw(ArgumentError(_alias_error_msg("surface", N)))
     end
 
     return integral(f, geometry, rule; kwargs...)
@@ -90,8 +92,7 @@ function volumeintegral(
         kwargs...
 )
     if (N = Meshes.paramdim(geometry)) != 3
-        throw(ArgumentError("Performing a volume integral on a geometry \
-                            with $N parametric dimensions not supported."))
+        throw(ArgumentError(_alias_error_msg("volume", N)))
     end
 
     return integral(f, geometry, rule; kwargs...)

src/integration_rules.jl

@@ -12,10 +12,9 @@ abstract type IntegrationRule end
     GaussKronrod(kwargs...)
 
 The h-adaptive Gauss-Kronrod quadrature rule implemented by
-[QuadGK.jl](https://github.com/JuliaMath/QuadGK.jl). All standard
-`QuadGK.quadgk` keyword arguments are supported. This rule works natively for one
-dimensional geometries; some two- and three-dimensional geometries are additionally
-supported using nested integral solvers with the specified `kwarg` settings.
+[QuadGK.jl](https://github.com/JuliaMath/QuadGK.jl) which can be used for any
+single-dimensional geometry. All standard `QuadGK.quadgk` keyword arguments are
+supported.
 """
 struct GaussKronrod <: IntegrationRule
     kwargs::Base.Pairs

src/utils.jl

@@ -1,13 +1,3 @@
-################################################################################
-#                           Misc. Internal Tools
-################################################################################
-
-# Common error message structure
-function _error_unsupported_combination(geometry, rule)
-    msg = "Integrating a $geometry using a $rule rule not supported."
-    throw(ArgumentError(msg))
-end
-
 ################################################################################
 #                           DifferentiationMethod
 ################################################################################

test/combinations.jl

@@ -63,7 +63,7 @@ This file includes tests for:
         elseif N == 2
             # surface
             aliases = Bool.((0, 1, 0))
-            rules = Bool.((1, 1, 1))
+            rules = Bool.((0, 1, 1))
             return SupportStatus(aliases..., rules..., autoenzyme)
         elseif N == 3
             # volume