From c4235b78a57bb9861e39dbbdf8ae0e8a2bdc6a1b Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 27 Nov 2024 17:01:41 -0500 Subject: [PATCH 01/14] Rely on automatic promotion instead of explicit conversion --- src/integral.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index 0fee71fb..c9f2164e 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -71,14 +71,14 @@ function _integral( nodes = Iterators.product(ntuple(Returns(xs), N)...) # Domain transformation: x [-1,1] ↦ t [0,1] - t(x) = FP(1 // 2) * x + FP(1 // 2) + t(x) = (1 // 2) * x + (1 // 2) function integrand((weights, nodes)) ts = t.(nodes) prod(weights) * f(geometry(ts...)) * differential(geometry, ts, diff_method) end - return FP(1 // (2^N)) .* sum(integrand, zip(weights, nodes)) + return (1 // (2^N)) .* sum(integrand, zip(weights, nodes)) end # HAdaptiveCubature From 54b8c390d3ee6fdc4c4b91abea192e59c82553b7 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 27 Nov 2024 17:27:15 -0500 Subject: [PATCH 02/14] Consolidate GaussKronrod code and issue depwarn in >1D --- src/integral.jl | 35 ++++++----------------------------- 1 file changed, 6 insertions(+), 29 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index c9f2164e..ec641cb4 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -46,9 +46,13 @@ function _integral( # Pass through to dim-specific workers in next section of this file N = Meshes.paramdim(geometry) if N == 1 - return _integral_gk_1d(f, geometry, rule; kwargs...) + integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method) + return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] elseif N == 2 - return _integral_gk_2d(f, geometry, rule; kwargs...) + Base.depwarn("Use `HAdaptiveCubature` instead of nested `GaussKronrod` rules.", :integral) + integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), diff_method) + ∫₁(v) = QuadGK.quadgk(u -> integrand(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") @@ -105,30 +109,3 @@ function _integral( # Reapply units return value .* integrandunits end - -################################################################################ -# Specialized GaussKronrod Methods -################################################################################ - -function _integral_gk_1d( - f, - geometry, - rule::GaussKronrod; - FP::Type{T} = Float64, - diff_method::DM = _default_method(geometry) -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - integrand(t) = f(geometry(t)) * differential(geometry, (t,), diff_method) - return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] -end - -function _integral_gk_2d( - f, - geometry2d, - rule::GaussKronrod; - FP::Type{T} = Float64, - diff_method::DM = _default_method(geometry2d) -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), diff_method) - ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1] - return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1] -end From 06aadad9d30169b006d7986d0bcdc6625d6de613 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 27 Nov 2024 18:12:56 -0500 Subject: [PATCH 03/14] Update support matrix docs --- docs/src/supportmatrix.md | 49 +++++++++++++++++++++------------------ 1 file changed, 26 insertions(+), 23 deletions(-) diff --git a/docs/src/supportmatrix.md b/docs/src/supportmatrix.md index 82fcd5f3..8aee0dbe 100644 --- a/docs/src/supportmatrix.md +++ b/docs/src/supportmatrix.md @@ -1,58 +1,61 @@ # Support Matrix -While this library aims to support all possible integration rules and [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) -geometry types, some combinations are ill-suited and some others are simply not yet -implemented. The following Support Matrix aims to capture the current development state of -all geometry/rule combinations. Entries with a green check mark are fully supported -and have passing unit tests that provide some confidence they produce accurate results. +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 an array of numerical integration rules and techniques. However, some +combinations of geomtry types and integration rules are ill-suited, and some others are simply +not yet yet implemented. The following Support Matrix captures the current state of support for +all geometry/rule combinations. Entries with a green check mark are fully supported and pass +unit tests designed to check for accuracy. In general, Gauss-Kronrod integration rules are recommended (and the default) for geometries -with one parametric dimension, e.g.: `Segment`, `BezierCurve`, and `Rope`. Gauss-Kronrod -rules can also be applied to some geometries with more dimensions by nesting multiple -integration solves, but this is inefficient. These Gauss-Kronrod rules are supported (but -not recommended) for surface-like geometries, but not for volume-like geometries. For -geometries with more than one parametric dimension, e.g. surfaces and volumes, H-Adaptive -Cubature integration rules are recommended (and the default). +with one parametric dimension, e.g.: `Segment`, `BezierCurve`, and `Rope`. or geometries with +more than one parametric dimension, e.g. surfaces and volumes, H-Adaptive Cubature rules are +recommended (and the default). + +While it is possible to apply nested Gauss-Kronrod rules to numerically integrate geometries +with more than one parametric dimension, this produces results that are strictly inferior to +using an equivalent H-Adaptive Cubature rule, so support for this usage is not recommended. | Symbol | Support Level | |--------|---------| -| ✅ | Supported, passes unit tests | +| ✅ | Supported | | 🎗️ | Planned to support in the future | +| ⚠️ | Deprecated | | 🛑 | Not supported | | `Meshes.Geometry` | Gauss-Legendre | Gauss-Kronrod | H-Adaptive Cubature | |----------|----------------|---------------|---------------------| -| `Ball` in `𝔼{2}` | ✅ | ✅ | ✅ | +| `Ball` in `𝔼{2}` | ✅ | ⚠️ | ✅ | | `Ball` in `𝔼{3}` | ✅ | 🛑 | ✅ | | `BezierCurve` | ✅ | ✅ | ✅ | | `Box` in `𝔼{1}` | ✅ | ✅ | ✅ | -| `Box` in `𝔼{2}` | ✅ | ✅ | ✅ | +| `Box` in `𝔼{2}` | ✅ | ⚠️ | ✅ | | `Box` in `𝔼{≥3}` | ✅ | 🛑 | ✅ | | `Circle` | ✅ | ✅ | ✅ | -| `Cone` | ✅ | ✅ | ✅ | -| `ConeSurface` | ✅ | ✅ | ✅ | +| `Cone` | ✅ | 🛑 | ✅ | +| `ConeSurface` | ✅ | ⚠️ | ✅ | | `Cylinder` | ✅ | 🛑 | ✅ | -| `CylinderSurface` | ✅ | ✅ | ✅ | -| `Disk` | ✅ | ✅ | ✅ | +| `CylinderSurface` | ✅ | ⚠️ | ✅ | +| `Disk` | ✅ | ⚠️ | ✅ | | `Ellipsoid` | ✅ | ✅ | ✅ | | `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) | -| `FrustumSurface` | ✅ | ✅ | ✅ | +| `FrustumSurface` | ✅ | ⚠️ | ✅ | | `Hexahedron` | ✅ | ✅ | ✅ | | `Line` | ✅ | ✅ | ✅ | -| `ParaboloidSurface` | ✅ | ✅ | ✅ | +| `ParaboloidSurface` | ✅ | ⚠️ | ✅ | | `ParametrizedCurve` | ✅ | ✅ | ✅ | | `Plane` | ✅ | ✅ | ✅ | | `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) | | `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) | -| `Quadrangle` | ✅ | ✅ | ✅ | +| `Quadrangle` | ✅ | ⚠️ | ✅ | | `Ray` | ✅ | ✅ | ✅ | | `Ring` | ✅ | ✅ | ✅ | | `Rope` | ✅ | ✅ | ✅ | | `Segment` | ✅ | ✅ | ✅ | | `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) | | `Sphere` in `𝔼{2}` | ✅ | ✅ | ✅ | -| `Sphere` in `𝔼{3}` | ✅ | ✅ | ✅ | +| `Sphere` in `𝔼{3}` | ✅ | ⚠️ | ✅ | | `Tetrahedron` in `𝔼{3}` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/40) | ✅ | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/40) | | `Triangle` | ✅ | ✅ | ✅ | -| `Torus` | ✅ | ✅ | ✅ | +| `Torus` | ✅ | ⚠️ | ✅ | | `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) | From 15a78f930c31f697c99ee2170569aa22cdf13eed Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 27 Nov 2024 19:03:55 -0500 Subject: [PATCH 04/14] Update function signature since not passing through --- src/integral.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index ec641cb4..d20e2008 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -41,8 +41,9 @@ function _integral( f, geometry, rule::GaussKronrod; - kwargs... -) + FP::Type{T} = Float64, + diff_method::DM = _default_method(geometry) +) where {DM <: DifferentiationMethod, T <: AbstractFloat} # Pass through to dim-specific workers in next section of this file N = Meshes.paramdim(geometry) if N == 1 From d69538dda08b0baf6ea98177010a18c7ee1057a3 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 27 Nov 2024 19:18:04 -0500 Subject: [PATCH 05/14] Update geometry naming --- src/integral.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integral.jl b/src/integral.jl index d20e2008..3301aef4 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -51,7 +51,7 @@ function _integral( return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] elseif N == 2 Base.depwarn("Use `HAdaptiveCubature` instead of nested `GaussKronrod` rules.", :integral) - integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v), diff_method) + integrand(u, v) = f(geometry(u, v)) * differential(geometry, (u, v), diff_method) ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1] return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1] else From 841e0d24978e35603a73e77e48aa50f9e0af9313 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 27 Nov 2024 19:35:19 -0500 Subject: [PATCH 06/14] Add comments and whitespace for clarity --- src/integral.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index 3301aef4..88d8a81b 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -44,16 +44,20 @@ function _integral( FP::Type{T} = Float64, diff_method::DM = _default_method(geometry) ) where {DM <: DifferentiationMethod, T <: AbstractFloat} - # Pass through to dim-specific workers in next section of this file N = Meshes.paramdim(geometry) + + # Implementation depends on number of parametric dimensions over which to integrate 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 nested `GaussKronrod` rules.", :integral) + + # Nested integration integrand(u, v) = f(geometry(u, v)) * differential(geometry, (u, v), diff_method) - ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1] - return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1] + ∫(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1] + return QuadGK.quadgk(∫, zero(FP), one(FP); rule.kwargs...)[1] else _error_unsupported_combination("geometry with more than two parametric dimensions", "GaussKronrod") From bd2f6509b5b45b547b61ba973615d4824e3d4eff Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 27 Nov 2024 19:36:14 -0500 Subject: [PATCH 07/14] Fix typo --- docs/src/supportmatrix.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/supportmatrix.md b/docs/src/supportmatrix.md index 8aee0dbe..7dc1b4ec 100644 --- a/docs/src/supportmatrix.md +++ b/docs/src/supportmatrix.md @@ -2,7 +2,7 @@ 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 an array of numerical integration rules and techniques. However, some -combinations of geomtry types and integration rules are ill-suited, and some others are simply +combinations of geometry types and integration rules are ill-suited, and some others are simply not yet yet implemented. The following Support Matrix captures the current state of support for all geometry/rule combinations. Entries with a green check mark are fully supported and pass unit tests designed to check for accuracy. From 57527c9cd7517d2f8a221b1c5d46e087ce5608f1 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 27 Nov 2024 20:42:12 -0500 Subject: [PATCH 08/14] Remove unused QuadGK dep --- test/Project.toml | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 81203128..e38225db 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,7 +4,6 @@ CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" -QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" @@ -17,7 +16,6 @@ CoordRefSystems = "0.12, 0.13, 0.14, 0.15, 0.16" ExplicitImports = "1.6.0" LinearAlgebra = "1" Meshes = "0.50, 0.51, 0.52" -QuadGK = "2.1.1" SpecialFunctions = "2" TestItemRunner = "1" TestItems = "1" From f2bf9ba9b3c3df176c0e6d1e6284c7b668205f7d Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 27 Nov 2024 20:55:25 -0500 Subject: [PATCH 09/14] Sanity check --- src/integral.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integral.jl b/src/integral.jl index 88d8a81b..551210ee 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -52,7 +52,7 @@ function _integral( return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1] elseif N == 2 # Issue deprecation warning - Base.depwarn("Use `HAdaptiveCubature` instead of nested `GaussKronrod` rules.", :integral) + #Base.depwarn("Use `HAdaptiveCubature` instead of nested `GaussKronrod` rules.", :integral) # Nested integration integrand(u, v) = f(geometry(u, v)) * differential(geometry, (u, v), diff_method) From 5eb6bfb02e515c770a2ea9d9f195dff1da7b7937 Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 27 Nov 2024 21:07:02 -0500 Subject: [PATCH 10/14] Sanity check part 2 --- src/integral.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index 551210ee..ed70f794 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -56,8 +56,8 @@ function _integral( # Nested integration integrand(u, v) = f(geometry(u, v)) * differential(geometry, (u, v), diff_method) - ∫(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1] - return QuadGK.quadgk(∫, zero(FP), one(FP); rule.kwargs...)[1] + ∫₁(v) = QuadGK.quadgk(u -> integrand(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") From ee5bdaea2d18926003f5ea16fa604c8bf54d493e Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 27 Nov 2024 21:27:36 -0500 Subject: [PATCH 11/14] Remove explicit FP conversion --- src/specializations/BezierCurve.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 7b0f44ab..c3540cbc 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -43,12 +43,12 @@ function integral( xs, ws = _gausslegendre(FP, rule.n) # Change of variables: x [-1,1] ↦ t [0,1] - t(x) = FP(1 // 2) * x + FP(1 // 2) + t(x) = (1 // 2) * x + (1 // 2) point(x) = curve(t(x), alg) integrand(x) = f(point(x)) * differential(curve, (t(x),), diff_method) # Integrate f along curve and apply domain-correction for [-1,1] ↦ [0, length] - return FP(1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs)) + return (1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs)) end function integral( From 4239cd23175dca6b274389be9c747254e185befe Mon Sep 17 00:00:00 2001 From: Mike Ingold Date: Wed, 27 Nov 2024 21:28:01 -0500 Subject: [PATCH 12/14] Add const for N's, and tweak formatting --- src/integral.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index ed70f794..e68030e4 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -33,7 +33,7 @@ function integral( end ################################################################################ -# Generalized (n-Dimensional) Worker Methods +# Integral Workers ################################################################################ # GaussKronrod @@ -44,15 +44,15 @@ function _integral( FP::Type{T} = Float64, diff_method::DM = _default_method(geometry) ) where {DM <: DifferentiationMethod, T <: AbstractFloat} - N = Meshes.paramdim(geometry) - # Implementation depends on number of parametric dimensions over which to integrate + const 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 nested `GaussKronrod` rules.", :integral) + Base.depwarn("Use `HAdaptiveCubature` instead of \ + `GaussKronrod` for surfaces.", :integral) # Nested integration integrand(u, v) = f(geometry(u, v)) * differential(geometry, (u, v), diff_method) @@ -72,7 +72,7 @@ function _integral( FP::Type{T} = Float64, diff_method::DM = _default_method(geometry) ) where {DM <: DifferentiationMethod, T <: AbstractFloat} - N = Meshes.paramdim(geometry) + const N = Meshes.paramdim(geometry) # Get Gauss-Legendre nodes and weights for a region [-1,1]^N xs, ws = _gausslegendre(FP, rule.n) @@ -98,7 +98,7 @@ function _integral( FP::Type{T} = Float64, diff_method::DM = _default_method(geometry) ) where {DM <: DifferentiationMethod, T <: AbstractFloat} - N = Meshes.paramdim(geometry) + const N = Meshes.paramdim(geometry) integrand(ts) = f(geometry(ts...)) * differential(geometry, ts, diff_method) From 5551192eed9371e7a397ee2af1a7e526b806318b Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Thu, 28 Nov 2024 08:54:20 -0500 Subject: [PATCH 13/14] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- docs/src/supportmatrix.md | 4 ++-- src/integral.jl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/src/supportmatrix.md b/docs/src/supportmatrix.md index 7dc1b4ec..4265c06a 100644 --- a/docs/src/supportmatrix.md +++ b/docs/src/supportmatrix.md @@ -3,12 +3,12 @@ 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 an array of numerical integration rules and techniques. However, some combinations of geometry types and integration rules are ill-suited, and some others are simply -not yet yet implemented. The following Support Matrix captures the current state of support for +not yet implemented. The following Support Matrix captures the current state of support for all geometry/rule combinations. Entries with a green check mark are fully supported and pass unit tests designed to check for accuracy. In general, Gauss-Kronrod integration rules are recommended (and the default) for geometries -with one parametric dimension, e.g.: `Segment`, `BezierCurve`, and `Rope`. or geometries with +with one parametric dimension, e.g.: `Segment`, `BezierCurve`, and `Rope`. For geometries with more than one parametric dimension, e.g. surfaces and volumes, H-Adaptive Cubature rules are recommended (and the default). diff --git a/src/integral.jl b/src/integral.jl index e68030e4..f6647e48 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -45,7 +45,7 @@ function _integral( diff_method::DM = _default_method(geometry) ) where {DM <: DifferentiationMethod, T <: AbstractFloat} # Implementation depends on number of parametric dimensions over which to integrate - const N = Meshes.paramdim(geometry) + 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] @@ -72,7 +72,7 @@ function _integral( FP::Type{T} = Float64, diff_method::DM = _default_method(geometry) ) where {DM <: DifferentiationMethod, T <: AbstractFloat} - const N = Meshes.paramdim(geometry) + N = Meshes.paramdim(geometry) # Get Gauss-Legendre nodes and weights for a region [-1,1]^N xs, ws = _gausslegendre(FP, rule.n) @@ -98,7 +98,7 @@ function _integral( FP::Type{T} = Float64, diff_method::DM = _default_method(geometry) ) where {DM <: DifferentiationMethod, T <: AbstractFloat} - const N = Meshes.paramdim(geometry) + N = Meshes.paramdim(geometry) integrand(ts) = f(geometry(ts...)) * differential(geometry, ts, diff_method) From f1e7d46fae4de3ccc88f003d5076afb59f2bd248 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Thu, 28 Nov 2024 08:54:49 -0500 Subject: [PATCH 14/14] Apply suggestion Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/integral.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/integral.jl b/src/integral.jl index f6647e48..de69db5d 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -55,8 +55,8 @@ function _integral( `GaussKronrod` for surfaces.", :integral) # Nested integration - integrand(u, v) = f(geometry(u, v)) * differential(geometry, (u, v), diff_method) - ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1] + 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",