Consolidate onto Jacobian-based integration methods

Project.toml

@@ -1,7 +1,7 @@
 name = "MeshIntegrals"
 uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
 authors = ["Mike Ingold <mike.ingold@gmail.com>"]
-version = "0.10.0"
+version = "0.11.0"
 
 [deps]
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"

README.md

@@ -30,7 +30,7 @@ Methods are tested to ensure compatibility with
 ### Integral
 | Geometry | Gauss-Legendre | Gauss-Kronrod | H-Adaptive Cubature |
 |----------|----------------|---------------|---------------------|
-| `Meshes.Box{Dim,T}` | :x: | :x: | :x: |
+| `Meshes.Box{Dim>3,T}` | :x: | :x: | :x: |
 
 ### Line Integral
 | Geometry | Gauss-Legendre | Gauss-Kronrod | H-Adaptive Cubature |
@@ -51,10 +51,10 @@ Methods are tested to ensure compatibility with
 | `Meshes.Box{2,T}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Meshes.CylinderSurface` | :x: | :white_check_mark: | :x: |
 | `Meshes.Disk` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
-| `Meshes.ParaboloidSurface` | :x: | :x: | :x: |
+| `Meshes.ParaboloidSurface` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `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.Torus` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
 | `Meshes.SimpleMesh` | :x: | :x: | :x: |
 
 ### Volume Integral
@@ -95,13 +95,15 @@ fr(p) = fr(p.coords...)
 
 # Plans and Work in Progress
 
-- Short term:
-    - Implement all methods in the support matrix above
+- TODO
     - Update Example Usage and benchmarks
     - Re-implement all tests for Unitful compatibility
-
-- Longer term:
-    - Once all methods are implemented, determine which can be consolidated with a more abstract/unioned case
     - Implement Aqua.jl tests
-    - Implement Documenter docs
+    - Improve README documentation in lieu of Documenter
     - Implement Monte Carlo integration methods
+    - Implement BenchmarkTools test suite to quantify performance
+    - Continue working to generalize and consolidate methods
+
+- Longer term plans
+    - Assess impact of upcoming Meshes.jl CRS refactor
+    - Once functionality is established and stable, evaluate transition to JuliaGeometry organization or direct absorption into Meshes.jl

src/integral_line.jl

@@ -34,6 +34,25 @@ end
 #                            Gauss-Legendre
 ################################################################################
 
+# Generalized method
+function _integral_1d(f, geometry, settings::GaussLegendre)
+    # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
+    xs, ws = gausslegendre(settings.n)
+
+    # Change of variables: x [-1,1] ↦ t [0,1]
+    t(x) = 0.5x + 0.5
+    point(x) = geometry(t(x))
+
+    function paramfactor(x)
+        J = jacobian(geometry,[t(x)])
+        return norm(J[1])
+    end
+
+    # Integrate f along the geometry and apply a domain-correction factor for [-1,1] ↦ [0, 1]
+    integrand((w,x)) = w * f(point(x)) * paramfactor(x)
+    return 0.5 * sum(integrand, zip(ws, xs))
+end
+
 """
     integral(f, curve::Meshes.BezierCurve, ::GaussLegendre; alg=Meshes.Horner())
 
@@ -73,15 +92,7 @@ function integral(
     # A Box{1,T} is definitionally embedded in 1D-space
     _validate_integrand(f,1,T)
 
-    # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
-    xs, ws = gausslegendre(settings.n)
-
-    # Change of variables: x [-1,1] ↦ t [0,1]
-    t(x) = 0.5x + 0.5
-    point(x) = line(t(x))
-
-    # Integrate f along the line and apply a domain-correction factor for [-1,1] ↦ [0, length]
-    return 0.5 * length(line) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs))
+    return _integral_1d(f, line, settings)
 end
 
 function integral(
@@ -93,16 +104,7 @@ function integral(
     # A Circle is definitionally embedded in 3D-space
     _validate_integrand(f,3,T)
 
-    # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
-    xs, ws = gausslegendre(settings.n)
-
-    # Change of variables: x [-1,1] ↦ t [0,1]
-    t(x) = 0.5x + 0.5
-    point(x) = circle(t(x))
-
-    # Integrate f along the circle's rim and apply a domain-correction
-    #   factor for [-1,1] ↦ [0, circumference]
-    return 0.5 * length(circle) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs))
+    return _integral_1d(f, circle, settings)
 end
 
 function integral(
@@ -135,15 +137,7 @@ function integral(
     # Validate the provided integrand function
     _validate_integrand(f,Dim,T)
 
-    # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
-    xs, ws = gausslegendre(settings.n)
-
-    # Change of variables: x [-1,1] ↦ t [0,1]
-    t(x) = 0.5x + 0.5
-    point(x) = segment(t(x))
-
-    # Integrate f along the line and apply a domain-correction factor for [-1,1] ↦ [0, length]
-    return 0.5 * length(segment) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs))
+    return _integral_1d(f, segment, settings)
 end
 
 function integral(
@@ -155,23 +149,25 @@ function integral(
     # A Sphere{2,T} is simply a circle in 2D-space
     _validate_integrand(f,2,T)
 
-    # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
-    xs, ws = gausslegendre(settings.n)
-
-    # Change of variables: x [-1,1] ↦ t [0,1]
-    t(x) = 0.5x + 0.5
-    point(x) = circle(t(x))
-
-    # Integrate f along the circle's rim and apply a domain-correction
-    #   factor for [-1,1] ↦ [0, circumference]
-    return 0.5 * length(circle) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs))
+    return _integral_1d(f, circle, settings)
 end
 
 
 ################################################################################
 #                                Gauss-Kronrod
 ################################################################################
 
+# Generalized method
+function _integral_1d(f, geometry, settings::GaussKronrod)
+    function paramfactor(t)
+        J = jacobian(geometry,[t])
+        return norm(J[1])
+    end
+
+    integrand(t) = f(geometry(t)) * paramfactor(t)
+    return QuadGK.quadgk(integrand, 0, 1; settings.kwargs...)[1]
+end
+
 """
     integral(f, curve::BezierCurve, ::GaussKronrod; alg=Horner(), kws...)
 
@@ -205,9 +201,7 @@ function integral(
     # A Box is definitionally embedded in 1D-space
     _validate_integrand(f,1,T)
 
-    len = length(line)
-    point(t) = line(t)
-    return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1]
+    return _integral_1d(f, line, settings)
 end
 
 function integral(
@@ -219,9 +213,7 @@ function integral(
     # A Circle is definitionally embedded in 3D-space
     _validate_integrand(f,3,T)
 
-    len = length(circle)
-    point(t) = circle(t)
-    return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1]
+    return _integral_1d(f, circle, settings)
 end
 
 function integral(
@@ -248,9 +240,7 @@ function integral(
     # Validate the provided integrand function
     _validate_integrand(f,Dim,T)
 
-    len = length(segment)
-    point(t) = segment(t)
-    return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1]
+    return _integral_1d(f, segment, settings)
 end
 
 function integral(
@@ -262,16 +252,25 @@ function integral(
     # A Sphere{2,T} is simply a circle in 2D-space
     _validate_integrand(f,2,T)
 
-    len = length(circle)
-    point(t) = circle(t)
-    return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1]
+    return _integral_1d(f, circle, settings)
 end
 
 
 ################################################################################
 #                                HCubature
 ################################################################################
 
+# Generalized method
+function _integral_1d(f, geometry, settings::HAdaptiveCubature)
+    function paramfactor(t)
+        J = jacobian(geometry,t)
+        return norm(J[1])
+    end
+
+    integrand(t) = f(geometry(t[1])) * paramfactor(t)
+    return HCubature.hcubature(integrand, [0], [1]; settings.kwargs...)[1]
+end
+
 """
     integral(f, curve::BezierCurve, ::HAdaptiveCubature; alg=Horner(), kws...)
 
@@ -305,9 +304,7 @@ function integral(
     # A Box is definitionally embedded in 1D-space
     _validate_integrand(f,1,T)
 
-    len = length(line)
-    point(t) = line(t)
-    return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1]
+    return _integral_1d(f, line, settings)
 end
 
 function integral(
@@ -319,9 +316,7 @@ function integral(
     # A Circle is definitionally embedded in 3D-space
     _validate_integrand(f,3,T)
 
-    len = length(circle)
-    point(t) = circle(t)
-    return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1]
+    return _integral_1d(f, circle, settings)
 end
 
 function integral(
@@ -352,9 +347,7 @@ function integral(
     # Validate the provided integrand function
     _validate_integrand(f,Dim,T)
 
-    len = length(segment)
-    point(t) = segment(t)
-    return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1]
+    return _integral_1d(f, segment, settings)
 end
 
 function integral(
@@ -366,7 +359,5 @@ function integral(
     # A Sphere{2,T} is simply a circle in 2D-space
     _validate_integrand(f,2,T)
 
-    len = length(circle)
-    point(t) = circle(t)
-    return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1]
+    return _integral_1d(f, circle, settings)
 end