Implement _ParametricGeometry (#138)

* Add a discrete specializations section

* Add Unitful and explicit julia compat

* Implement _ParametricGeometry and use it for BezierCurve and Tetrahedron

* Update Tetrahedron tests

* Update tests

* Update a remnant usage of zeros

* Fix Unitful integrand handling in Triangle HAdaptiveCubature method

* Add missing reapplication of units

* Remove extended tag from Tetrahedron tests

* Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Add tests for _parametric bounds checks

* Fix typo

* Fix typo

* Add checks for t1 and t2 args of _parametric Tetrahedron

* Update Tetrahedron status

* Apply format suggestion

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Apply suggestions from code review

Co-authored-by: Joshua Lampert <[email protected]>

* Update src/specializations/BezierCurve.jl

Co-authored-by: Joshua Lampert <[email protected]>

* Update src/specializations/Tetrahedron.jl

Co-authored-by: Joshua Lampert <[email protected]>

* Apply suggestions from code review

Co-authored-by: Joshua Lampert <[email protected]>

---------

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: Joshua Lampert <[email protected]>

Author: 3 people

benchmark/Project.toml

@@ -2,8 +2,11 @@
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 BenchmarkTools = "1.5"
 LinearAlgebra = "1"
 Meshes = "0.50, 0.51, 0.52"
+Unitful = "1.19"
+julia = "1.9"

benchmark/benchmarks.jl

@@ -2,6 +2,7 @@ using BenchmarkTools
 using LinearAlgebra
 using Meshes
 using MeshIntegrals
+using Unitful
 
 const SUITE = BenchmarkGroup()
 
@@ -19,10 +20,8 @@ rules = (
     (name = "HAdaptiveCubature", rule = HAdaptiveCubature(), N = 500)
 )
 geometries = (
-    (name = "Meshes.BezierCurve",
-        item = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)])),
-    (name = "Meshes.Segment", item = Segment(Point(0, 0, 0), Point(1, 1, 1))),
-    (name = "Meshes.Sphere", item = Sphere(Point(0, 0, 0), 1.0))
+    (name = "Segment", item = Segment(Point(0, 0, 0), Point(1, 1, 1))),
+    (name = "Sphere", item = Sphere(Point(0, 0, 0), 1.0))
 )
 
 SUITE["Integrals"] = let s = BenchmarkGroup()
@@ -33,6 +32,43 @@ SUITE["Integrals"] = let s = BenchmarkGroup()
     s
 end
 
+############################################################################################
+#                                    Specializations
+############################################################################################
+
+spec = (
+    f = p -> norm(to(p)),
+    f_exp = p::Point -> exp(-norm(to(p))^2 / u"m^2"),
+    g = (
+        bezier = BezierCurve([Point(t, sin(t), 0) for t in range(-pi, pi, length = 361)]),
+        line = Line(Point(0, 0, 0), Point(1, 1, 1)),
+        plane = Plane(Point(0, 0, 0), Vec(0, 0, 1)),
+        ray = Ray(Point(0, 0, 0), Vec(0, 0, 1)),
+        triangle = Triangle(Point(1, 0, 0), Point(0, 1, 0), Point(0, 0, 1)),
+        tetrahedron = let
+            a = Point(0, 3, 0)
+            b = Point(-7, 0, 0)
+            c = Point(8, 0, 0)
+            ẑ = Vec(0, 0, 1)
+            Tetrahedron(a, b, c, a + ẑ)
+        end
+    ),
+    rule_gl = GaussLegendre(100),
+    rule_gk = GaussKronrod(),
+    rule_hc = HAdaptiveCubature()
+)
+
+SUITE["Specializations/Scalar GaussLegendre"] = let s = BenchmarkGroup()
+    s["BezierCurve"] = @benchmarkable integral($spec.f, $spec.g.bezier, $spec.rule_gl)
+    s["Line"] = @benchmarkable integral($spec.f_exp, $spec.g.line, $spec.rule_gl)
+    s["Plane"] = @benchmarkable integral($spec.f_exp, $spec.g.plane, $spec.rule_gl)
+    s["Ray"] = @benchmarkable integral($spec.f_exp, $spec.g.ray, $spec.rule_gl)
+    s["Triangle"] = @benchmarkable integral($spec.f, $spec.g.triangle, $spec.rule_gl)
+    # s["Tetrahedron"] = @benchmarkable integral($spec.f, $spec.g.tetrahedron, $spec.rule)
+    # TODO re-enable Tetrahedron-GaussLegendre test when supported by main branch
+    s
+end
+
 ############################################################################################
 #                                      Differentials
 ############################################################################################

docs/src/supportmatrix.md

@@ -55,7 +55,7 @@ using an equivalent H-Adaptive Cubature rule, so support for this usage is not r
 | `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}` | ✅ | ⚠️ | ✅ |
-| `Tetrahedron` in `𝔼{3}` | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/40) | ✅ | [🎗️](https://github.com/JuliaGeometry/MeshIntegrals.jl/issues/40) |
+| `Tetrahedron` | ✅ | ⚠️ | ✅ |
 | `Triangle` | ✅ | ✅ | ✅ |
 | `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) |

src/MeshIntegrals.jl

@@ -24,6 +24,7 @@ include("integral_aliases.jl")
 export lineintegral, surfaceintegral, volumeintegral
 
 # Integration methods specialized for particular geometries
+include("specializations/_ParametricGeometry.jl")
 include("specializations/BezierCurve.jl")
 include("specializations/ConeSurface.jl")
 include("specializations/CylinderSurface.jl")

src/integral.jl

@@ -104,7 +104,7 @@ function _integral(
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
-    testpoint_parametriccoord = zeros(FP, N)
+    testpoint_parametriccoord = _zeros(FP, N)
     integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
     # Create a wrapper that returns only the value component in those units
     uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts))

src/specializations/BezierCurve.jl

@@ -34,93 +34,22 @@ steep performance cost, especially for curves with a large number of control poi
 function integral(
         f,
         curve::Meshes.BezierCurve,
-        rule::GaussLegendre;
-        diff_method::DM = _default_method(curve),
-        FP::Type{T} = Float64,
-        alg::Meshes.BezierEvalMethod = Meshes.Horner()
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
-    xs, ws = _gausslegendre(FP, rule.n)
-
-    # Change of variables: x [-1,1] ↦ t [0,1]
-    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 (1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs))
-end
-
-function integral(
-        f,
-        curve::Meshes.BezierCurve,
-        rule::GaussKronrod;
-        diff_method::DM = _default_method(curve),
-        FP::Type{T} = Float64,
-        alg::Meshes.BezierEvalMethod = Meshes.Horner()
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    point(t) = curve(t, alg)
-    integrand(t) = f(point(t)) * differential(curve, (t,), diff_method)
-    return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
-end
-
-function integral(
-        f,
-        curve::Meshes.BezierCurve,
-        rule::HAdaptiveCubature;
-        diff_method::DM = _default_method(curve),
-        FP::Type{T} = Float64,
-        alg::Meshes.BezierEvalMethod = Meshes.Horner()
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    point(t) = curve(t, alg)
-    integrand(ts) = f(point(only(ts))) * differential(curve, ts, diff_method)
-
-    # HCubature doesn't support functions that output Unitful Quantity types
-    # Establish the units that are output by f
-    testpoint_parametriccoord = zeros(FP, 1)
-    integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
-    # Create a wrapper that returns only the value component in those units
-    uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts))
-    # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, _zeros(FP, 1), _ones(FP, 1); rule.kwargs...)[1]
-
-    # Reapply units
-    return value .* integrandunits
+        rule::IntegrationRule;
+        alg::Meshes.BezierEvalMethod = Meshes.Horner(),
+        kwargs...
+)
+    # Generate a _ParametricGeometry whose parametric function auto-applies the alg kwarg
+    paramfunction(t) = _parametric(curve, t, alg)
+    param_curve = _ParametricGeometry(paramfunction, Meshes.paramdim(curve))
+
+    # Integrate the _ParametricGeometry using the standard methods
+    return _integral(f, param_curve, rule; kwargs...)
 end
 
 ################################################################################
-#                               jacobian
+#                              Parametric
 ################################################################################
 
-function jacobian(
-        bz::Meshes.BezierCurve,
-        ts::V,
-        diff_method::Analytical
-) where {V <: Union{AbstractVector, Tuple}}
-    t = only(ts)
-    # Parameter t restricted to domain [0,1] by definition
-    if t < 0 || t > 1
-        throw(DomainError(t, "b(t) is not defined for t outside [0, 1]."))
-    end
-
-    # Aliases
-    P = bz.controls
-    N = Meshes.degree(bz)
-
-    # Ensure that this implementation is tractible: limited by ability to calculate
-    #   binomial(N, N/2) without overflow. It's possible to extend this range by
-    #   converting N to a BigInt, but this results in always returning BigFloat's.
-    N <= 1028 || error("This algorithm overflows for curves with ⪆1000 control points.")
-
-    # Generator for Bernstein polynomial functions
-    B(i, n) = t -> binomial(Int128(n), i) * t^i * (1 - t)^(n - i)
-
-    # Derivative = N Σ_{i=0}^{N-1} sigma(i)
-    #   P indices adjusted for Julia 1-based array indexing
-    sigma(i) = B(i, N - 1)(t) .* (P[(i + 1) + 1] - P[(i) + 1])
-    derivative = N .* sum(sigma, 0:(N - 1))
-
-    return (derivative,)
+function _parametric(curve::Meshes.BezierCurve, t, alg::Meshes.BezierEvalMethod)
+    return curve(t, alg)
 end
-
-_has_analytical(::Type{T}) where {T <: Meshes.BezierCurve} = true

src/specializations/Tetrahedron.jl

@@ -12,44 +12,32 @@
 function integral(
         f,
         tetrahedron::Meshes.Tetrahedron,
-        rule::GaussLegendre;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _error_unsupported_combination("Tetrahedron", "GaussLegendre")
-end
-
-function integral(
-        f,
-        tetrahedron::Meshes.Tetrahedron,
-        rule::GaussKronrod;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Tetrahedron, diff_method)
-
-    o = zero(FP)
-    ∫uvw(u, v, w) = f(tetrahedron(u, v, w))
-    ∫vw(v, w) = QuadGK.quadgk(u -> ∫uvw(u, v, w), o, FP(1 - v - w); rule.kwargs...)[1]
-    ∫w(w) = QuadGK.quadgk(v -> ∫vw(v, w), o, FP(1 - w); rule.kwargs...)[1]
-    ∫ = QuadGK.quadgk(∫w, o, one(FP); rule.kwargs...)[1]
-
-    # Apply barycentric domain correction (volume: 1/6 → actual)
-    return 6 * Meshes.volume(tetrahedron) * ∫
-end
+        rule::IntegrationRule;
+        kwargs...
+)
+    # Generate a _ParametricGeometry whose parametric function domain spans [0,1]^3
+    paramfunction(t1, t2, t3) = _parametric(tetrahedron, t1, t2, t3)
+    tetra = _ParametricGeometry(paramfunction, Meshes.paramdim(tetrahedron))
 
-function integral(
-        f,
-        tetrahedron::Meshes.Tetrahedron,
-        rule::HAdaptiveCubature;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _error_unsupported_combination("Tetrahedron", "HAdaptiveCubature")
+    # Integrate the _ParametricGeometry using the standard methods
+    return _integral(f, tetra, rule; kwargs...)
 end
 
 ################################################################################
-#                               jacobian
+#                               Parametric
 ################################################################################
 
-_has_analytical(::Type{T}) where {T <: Meshes.Tetrahedron} = true
+function _parametric(tetrahedron::Meshes.Tetrahedron, t1, t2, t3)
+    if (t3 < 0) || (t3 > 1)
+        msg = "tetrahedron(t1, t2, t3) is not defined for t3 outside [0, 1]."
+        throw(DomainError(t3, msg))
+    end
+
+    # Take a triangular cross-section at t3
+    a = tetrahedron(t3, 0, 0)
+    b = tetrahedron(0, t3, 0)
+    c = tetrahedron(0, 0, t3)
+    cross_section = Meshes.Triangle(a, b, c)
+
+    return _parametric(cross_section, t1, t2)
+end

src/specializations/Triangle.jl

@@ -87,14 +87,40 @@ function integral(
         v = R * (1 - b / (a + b))
         return f(triangle(u, v)) * R / (a + b)^2
     end
-    ∫ = HCubature.hcubature(integrand, _zeros(FP, 2), (FP(1), FP(π / 2)), rule.kwargs...)[1]
+
+    # HCubature doesn't support functions that output Unitful Quantity types
+    # Establish the units that are output by f
+    testpoint_parametriccoord = _zeros(T, 2)
+    integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
+    # Create a wrapper that returns only the value component in those units
+    uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts))
+    # Integrate only the unitless values
+    ∫ = HCubature.hcubature(uintegrand, _zeros(FP, 2), (T(1), T(π / 2)), rule.kwargs...)[1]
 
     # Apply a linear domain-correction factor 0.5 ↦ area(triangle)
-    return 2 * Meshes.area(triangle) .* ∫
+    return 2 * Meshes.area(triangle) * integrandunits .* ∫
 end
 
 ################################################################################
 #                               jacobian
 ################################################################################
 
 _has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true
+
+################################################################################
+#                              Parametric
+################################################################################
+
+function _parametric(triangle::Meshes.Triangle, t1, t2)
+    if (t1 < 0 || t1 > 1) || (t2 < 0 || t2 > 1)
+        msg = "triangle(t1, t2) is not defined for (t1, t2) outside [0, 1]^2."
+        throw(DomainError((t1, t2), msg))
+    end
+
+    # Form a line segment between sides
+    a = triangle(0, t2)
+    b = triangle(t2, 0)
+    segment = Meshes.Segment(a, b)
+
+    return segment(t1)
+end

src/specializations/_ParametricGeometry.jl

@@ -0,0 +1,63 @@
+"""
+    _ParametricGeometry <: Meshes.Primitive <: Meshes.Geometry
+
+_ParametricGeometry is used internally in MeshIntegrals.jl to behave like a generic wrapper
+for geometries with custom parametric functions. This type is used for transforming other
+geometries to enable integration over the standard rectangular `[0,1]^n` domain.
+
+Meshes.jl adopted a `ParametrizedCurve` type that performs a similar role as of `v0.51.20`,
+but only supports geometries with one parametric dimension. Support is additionally planned
+for more types that span surfaces and volumes, at which time this custom type will probably
+no longer be required.
+
+# Fields
+- `fun` - anything callable representing a parametric function: (ts...) -> Meshes.Point
+
+# Type Structure
+- `M <: Meshes.Manifold` - same usage as in `Meshes.Geometry{M, C}`
+- `C <: CoordRefSystems.CRS` - same usage as in `Meshes.Geometry{M, C}`
+- `F` - type of the callable parametric function
+- `Dim` - number of parametric dimensions
+"""
+struct _ParametricGeometry{M <: Meshes.Manifold, C <: CRS, F, Dim} <:
+       Meshes.Primitive{M, C}
+    fun::F
+
+    function _ParametricGeometry{M, C}(
+            fun::F,
+            Dim::Int64
+    ) where {M <: Meshes.Manifold, C <: CRS, F}
+        return new{M, C, F, Dim}(fun)
+    end
+end
+
+"""
+    _ParametricGeometry(fun, dims)
+
+Construct a `_ParametricGeometry` using a provided parametric function `fun` for a geometry
+with `dims` parametric dimensions.
+
+# Arguments
+- `fun` - anything callable representing a parametric function mapping `(ts...) -> Meshes.Point`
+- `dims::Int64` - number of parametric dimensions, i.e. `length(ts)`
+"""
+function _ParametricGeometry(
+        fun,
+        Dim::Int64
+)
+    p = fun(_zeros(Dim)...)
+    return _ParametricGeometry{Meshes.manifold(p), Meshes.crs(p)}(fun, Dim)
+end
+
+# Allow a _ParametricGeometry to be called like a Geometry
+(g::_ParametricGeometry)(ts...) = g.fun(ts...)
+
+Meshes.paramdim(::_ParametricGeometry{M, C, F, Dim}) where {M, C, F, Dim} = Dim
+
+"""
+    _parametric(geometry::G, ts...) where {G <: Meshes.Geometry}
+
+Used in MeshIntegrals.jl for defining parametric functions that transform non-standard
+geometries into a form that can be integrated over the standard rectangular [0,1]^n limits.
+"""
+function _parametric end