Skip to content

Commit cba39f0

Browse files
mikeingoldJoshuaLampert
mikeingold
and
JoshuaLampert
authored
Deprecate applying Gauss-Kronrod rules to surfaces (#136)
* Rely on automatic promotion instead of explicit conversion * Consolidate GaussKronrod code and issue depwarn in >1D * Update support matrix docs * Update function signature since not passing through * Update geometry naming * Add comments and whitespace for clarity * Fix typo * Remove unused QuadGK dep * Sanity check * Sanity check part 2 * Remove explicit FP conversion * Add const for N's, and tweak formatting * Apply suggestions from code review Co-authored-by: Joshua Lampert <[email protected]> * Apply suggestion Co-authored-by: Joshua Lampert <[email protected]> --------- Co-authored-by: Joshua Lampert <[email protected]>
1 parent 0cd335d commit cba39f0

File tree

4 files changed

+45
-62
lines changed

4 files changed

+45
-62
lines changed

‎docs/src/supportmatrix.md

Lines changed: 26 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,58 +1,61 @@
11
# Support Matrix
22

3-
While this library aims to support all possible integration rules and [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl)
4-
geometry types, some combinations are ill-suited and some others are simply not yet
5-
implemented. The following Support Matrix aims to capture the current development state of
6-
all geometry/rule combinations. Entries with a green check mark are fully supported
7-
and have passing unit tests that provide some confidence they produce accurate results.
3+
This library aims to enable users to calculate the value of integrals over all [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl)
4+
geometry types using an array of numerical integration rules and techniques. However, some
5+
combinations of geometry types and integration rules are ill-suited, and some others are simply
6+
not yet implemented. The following Support Matrix captures the current state of support for
7+
all geometry/rule combinations. Entries with a green check mark are fully supported and pass
8+
unit tests designed to check for accuracy.
89

910
In general, Gauss-Kronrod integration rules are recommended (and the default) for geometries
10-
with one parametric dimension, e.g.: `Segment`, `BezierCurve`, and `Rope`. Gauss-Kronrod
11-
rules can also be applied to some geometries with more dimensions by nesting multiple
12-
integration solves, but this is inefficient. These Gauss-Kronrod rules are supported (but
13-
not recommended) for surface-like geometries, but not for volume-like geometries. For
14-
geometries with more than one parametric dimension, e.g. surfaces and volumes, H-Adaptive
15-
Cubature integration rules are recommended (and the default).
11+
with one parametric dimension, e.g.: `Segment`, `BezierCurve`, and `Rope`. For geometries with
12+
more than one parametric dimension, e.g. surfaces and volumes, H-Adaptive Cubature rules are
13+
recommended (and the default).
14+
15+
While it is possible to apply nested Gauss-Kronrod rules to numerically integrate geometries
16+
with more than one parametric dimension, this produces results that are strictly inferior to
17+
using an equivalent H-Adaptive Cubature rule, so support for this usage is not recommended.
1618

1719
| Symbol | Support Level |
1820
|--------|---------|
19-
| ✅ | Supported, passes unit tests |
21+
| ✅ | Supported |
2022
| 🎗️ | Planned to support in the future |
23+
| ⚠️ | Deprecated |
2124
| 🛑 | Not supported |
2225

2326
| `Meshes.Geometry` | Gauss-Legendre | Gauss-Kronrod | H-Adaptive Cubature |
2427
|----------|----------------|---------------|---------------------|
25-
| `Ball` in `𝔼{2}` | ✅ | ✅ | ✅ |
28+
| `Ball` in `𝔼{2}` | ✅ | ⚠️ | ✅ |
2629
| `Ball` in `𝔼{3}` | ✅ | 🛑 | ✅ |
2730
| `BezierCurve` | ✅ | ✅ | ✅ |
2831
| `Box` in `𝔼{1}` | ✅ | ✅ | ✅ |
29-
| `Box` in `𝔼{2}` | ✅ | ✅ | ✅ |
32+
| `Box` in `𝔼{2}` | ✅ | ⚠️ | ✅ |
3033
| `Box` in `𝔼{≥3}` | ✅ | 🛑 | ✅ |
3134
| `Circle` | ✅ | ✅ | ✅ |
32-
| `Cone` | ✅ | ✅ | ✅ |
33-
| `ConeSurface` | ✅ | ✅ | ✅ |
35+
| `Cone` | ✅ | 🛑 | ✅ |
36+
| `ConeSurface` | ✅ | ⚠️ | ✅ |
3437
| `Cylinder` | ✅ | 🛑 | ✅ |
35-
| `CylinderSurface` | ✅ | ✅ | ✅ |
36-
| `Disk` | ✅ | ✅ | ✅ |
38+
| `CylinderSurface` | ✅ | ⚠️ | ✅ |
39+
| `Disk` | ✅ | ⚠️ | ✅ |
3740
| `Ellipsoid` | ✅ | ✅ | ✅ |
3841
| `Frustum` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) |
39-
| `FrustumSurface` | ✅ | ✅ | ✅ |
42+
| `FrustumSurface` | ✅ | ⚠️ | ✅ |
4043
| `Hexahedron` | ✅ | ✅ | ✅ |
4144
| `Line` | ✅ | ✅ | ✅ |
42-
| `ParaboloidSurface` | ✅ | ✅ | ✅ |
45+
| `ParaboloidSurface` | ✅ | ⚠️ | ✅ |
4346
| `ParametrizedCurve` | ✅ | ✅ | ✅ |
4447
| `Plane` | ✅ | ✅ | ✅ |
4548
| `Polyarea` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) |
4649
| `Pyramid` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) |
47-
| `Quadrangle` | ✅ | ✅ | ✅ |
50+
| `Quadrangle` | ✅ | ⚠️ | ✅ |
4851
| `Ray` | ✅ | ✅ | ✅ |
4952
| `Ring` | ✅ | ✅ | ✅ |
5053
| `Rope` | ✅ | ✅ | ✅ |
5154
| `Segment` | ✅ | ✅ | ✅ |
5255
| `SimpleMesh` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/27) |
5356
| `Sphere` in `𝔼{2}` | ✅ | ✅ | ✅ |
54-
| `Sphere` in `𝔼{3}` | ✅ | ✅ | ✅ |
57+
| `Sphere` in `𝔼{3}` | ✅ | ⚠️ | ✅ |
5558
| `Tetrahedron` in `𝔼{3}` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/40) | ✅ | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/40) |
5659
| `Triangle` | ✅ | ✅ | ✅ |
57-
| `Torus` | ✅ | ✅ | ✅ |
60+
| `Torus` | ✅ | ⚠️ | ✅ |
5861
| `Wedge` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/28) |

‎src/integral.jl

Lines changed: 17 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -33,22 +33,31 @@ function integral(
3333
end
3434

3535
################################################################################
36-
# Generalized (n-Dimensional) Worker Methods
36+
# Integral Workers
3737
################################################################################
3838

3939
# GaussKronrod
4040
function _integral(
4141
f,
4242
geometry,
4343
rule::GaussKronrod;
44-
kwargs...
45-
)
46-
# Pass through to dim-specific workers in next section of this file
44+
FP::Type{T} = Float64,
45+
diff_method::DM = _default_method(geometry)
46+
) where {DM <: DifferentiationMethod, T <: AbstractFloat}
47+
# Implementation depends on number of parametric dimensions over which to integrate
4748
N = Meshes.paramdim(geometry)
4849
if N == 1
49-
return _integral_gk_1d(f, geometry, rule; kwargs...)
50+
integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method)
51+
return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
5052
elseif N == 2
51-
return _integral_gk_2d(f, geometry, rule; kwargs...)
53+
# Issue deprecation warning
54+
Base.depwarn("Use `HAdaptiveCubature` instead of \
55+
`GaussKronrod` for surfaces.", :integral)
56+
57+
# Nested integration
58+
integrand2D(u, v) = f(geometry(u, v)) * differential(geometry, (u, v), diff_method)
59+
∫₁(v) = QuadGK.quadgk(u -> integrand2D(u, v), zero(FP), one(FP); rule.kwargs...)[1]
60+
return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1]
5261
else
5362
_error_unsupported_combination("geometry with more than two parametric dimensions",
5463
"GaussKronrod")
@@ -71,14 +80,14 @@ function _integral(
7180
nodes = Iterators.product(ntuple(Returns(xs), N)...)
7281

7382
# Domain transformation: x [-1,1] ↦ t [0,1]
74-
t(x) = FP(1 // 2) * x + FP(1 // 2)
83+
t(x) = (1 // 2) * x + (1 // 2)
7584

7685
function integrand((weights, nodes))
7786
ts = t.(nodes)
7887
prod(weights) * f(geometry(ts...)) * differential(geometry, ts, diff_method)
7988
end
8089

81-
return FP(1 // (2^N)) .* sum(integrand, zip(weights, nodes))
90+
return (1 // (2^N)) .* sum(integrand, zip(weights, nodes))
8291
end
8392

8493
# HAdaptiveCubature
@@ -105,30 +114,3 @@ function _integral(
105114
# Reapply units
106115
return value .* integrandunits
107116
end
108-
109-
################################################################################
110-
# Specialized GaussKronrod Methods
111-
################################################################################
112-
113-
function _integral_gk_1d(
114-
f,
115-
geometry,
116-
rule::GaussKronrod;
117-
FP::Type{T} = Float64,
118-
diff_method::DM = _default_method(geometry)
119-
) where {DM <: DifferentiationMethod, T <: AbstractFloat}
120-
integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method)
121-
return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
122-
end
123-
124-
function _integral_gk_2d(
125-
f,
126-
geometry2d,
127-
rule::GaussKronrod;
128-
FP::Type{T} = Float64,
129-
diff_method::DM = _default_method(geometry2d)
130-
) where {DM <: DifferentiationMethod, T <: AbstractFloat}
131-
integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), diff_method)
132-
∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1]
133-
return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1]
134-
end

‎src/specializations/BezierCurve.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -43,12 +43,12 @@ function integral(
4343
xs, ws = _gausslegendre(FP, rule.n)
4444

4545
# Change of variables: x [-1,1] ↦ t [0,1]
46-
t(x) = FP(1 // 2) * x + FP(1 // 2)
46+
t(x) = (1 // 2) * x + (1 // 2)
4747
point(x) = curve(t(x), alg)
4848
integrand(x) = f(point(x)) * differential(curve, (t(x),), diff_method)
4949

5050
# Integrate f along curve and apply domain-correction for [-1,1] ↦ [0, length]
51-
return FP(1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs))
51+
return (1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs))
5252
end
5353

5454
function integral(

‎test/Project.toml

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,6 @@ CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
44
ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
55
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
66
Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
7-
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
87
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
98
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
109
TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
@@ -17,7 +16,6 @@ CoordRefSystems = "0.12, 0.13, 0.14, 0.15, 0.16"
1716
ExplicitImports = "1.6.0"
1817
LinearAlgebra = "1"
1918
Meshes = "0.50, 0.51, 0.52"
20-
QuadGK = "2.1.1"
2119
SpecialFunctions = "2"
2220
TestItemRunner = "1"
2321
TestItems = "1"

0 commit comments

Comments
 (0)