Implement more _ParametricGeometry transforms

.typos.toml

@@ -0,0 +1,2 @@
+[default.extend-words]
+parametrized = "parametrized"

benchmark/benchmarks.jl

@@ -15,9 +15,9 @@ integrands = (
     (name = "Vector", f = p -> fill(norm(to(p)), 3))
 )
 rules = (
-    (name = "GaussLegendre", rule = GaussLegendre(100), N = 100),
-    (name = "GaussKronrod", rule = GaussKronrod(), N = 100),
-    (name = "HAdaptiveCubature", rule = HAdaptiveCubature(), N = 500)
+    (name = "GaussLegendre", rule = GaussLegendre(100)),
+    (name = "GaussKronrod", rule = GaussKronrod()),
+    (name = "HAdaptiveCubature", rule = HAdaptiveCubature())
 )
 geometries = (
     (name = "Segment", item = Segment(Point(0, 0, 0), Point(1, 1, 1))),
@@ -26,7 +26,8 @@ geometries = (
 
 SUITE["Integrals"] = let s = BenchmarkGroup()
     for (int, rule, geometry) in Iterators.product(integrands, rules, geometries)
-        n1, n2, N = geometry.name, "$(int.name) $(rule.name)", rule.N
+        n1 = geometry.name
+        n2 = "$(int.name) $(rule.name)"
         s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule)
     end
     s
@@ -64,8 +65,7 @@ SUITE["Specializations/Scalar GaussLegendre"] = let s = BenchmarkGroup()
     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["Tetrahedron"] = @benchmarkable integral($spec.f, $spec.g.tetrahedron, $spec.rule_gl)
     s
 end
 

src/specializations/BezierCurve.jl

@@ -39,8 +39,7 @@ function integral(
         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))
+    param_curve = _ParametricGeometry(_parametric(curve, alg), Meshes.paramdim(curve))
 
     # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_curve, rule; kwargs...)
@@ -50,6 +49,7 @@ end
 #                              Parametric
 ################################################################################
 
-function _parametric(curve::Meshes.BezierCurve, t, alg::Meshes.BezierEvalMethod)
-    return curve(t, alg)
+# Wrap (::BezierCurve)(t, ::BezierEvalMethod) into f(t) by embedding second argument
+function _parametric(curve::Meshes.BezierCurve, alg::Meshes.BezierEvalMethod)
+    return t -> curve(t, alg)
 end

src/specializations/Line.jl

@@ -10,81 +10,26 @@
 function integral(
         f,
         line::Meshes.Line,
-        rule::GaussLegendre;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Line, diff_method)
-
-    # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
-    xs, ws = _gausslegendre(FP, rule.n)
-
-    # Normalize the Line s.t. line(t) is distance t from origin point
-    line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a))
-
-    # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞)
-    t(x) = x / (1 - x^2)
-    t′(x) = (1 + x^2) / (1 - x^2)^2
-
-    # Integrate f along the Line
-    differential(line, x) = t′(x) * _units(line(0))
-    integrand(x) = f(line(t(x))) * differential(line, x)
-    return sum(w .* integrand(x) for (w, x) in zip(ws, xs))
-end
-
-function integral(
-        f,
-        line::Meshes.Line,
-        rule::GaussKronrod;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Line, diff_method)
-
-    # Normalize the Line s.t. line(t) is distance t from origin point
-    line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a))
-
-    # Integrate f along the Line
-    domainunits = _units(line(0))
-    integrand(t) = f(line(t)) * domainunits
-    return QuadGK.quadgk(integrand, FP(-Inf), FP(Inf); rule.kwargs...)[1]
-end
-
-function integral(
-        f,
-        line::Meshes.Line,
-        rule::HAdaptiveCubature;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Line, diff_method)
-
-    # Normalize the Line s.t. line(t) is distance t from origin point
-    line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a))
-
-    # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞)
-    t(x) = x / (1 - x^2)
-    t′(x) = (1 + x^2) / (1 - x^2)^2
-
-    # Integrate f along the Line
-    differential(line, x) = t′(x) * _units(line(0))
-    integrand(xs) = f(line(t(xs[1]))) * differential(line, xs[1])
-
-    # HCubature doesn't support functions that output Unitful Quantity types
-    # Establish the units that are output by f
-    testpoint_parametriccoord = (FP(0.5),)
-    integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
-    # Create a wrapper that returns only the value component in those units
-    uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
-    # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, (-one(FP),), (one(FP),); rule.kwargs...)[1]
-
-    # Reapply units
-    return value .* integrandunits
+        rule::IntegrationRule;
+        kwargs...
+)
+    # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]
+    param_line = _ParametricGeometry(_parametric(line), Meshes.paramdim(line))
+
+    # Integrate the _ParametricGeometry using the standard methods
+    return _integral(f, param_line, rule; kwargs...)
 end
 
 ################################################################################
-#                               jacobian
+#                              Parametric
 ################################################################################
 
-_has_analytical(::Type{T}) where {T <: Meshes.Line} = true
+# Map argument domain from [0, 1] to (-∞, ∞) for (::Line)(t)
+function _parametric(line::Meshes.Line)
+    # [-1, 1] ↦ (-∞, ∞)
+    f1(t) = t / (1 - t^2)
+    # [0, 1] ↦ [-1, 1]
+    f2(t) = 2t - 1
+    # Compose the two transforms
+    return t -> line((f1 ∘ f2)(t))
+end

src/specializations/Plane.jl

@@ -11,96 +11,27 @@
 function integral(
         f,
         plane::Meshes.Plane,
-        rule::GaussLegendre;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Plane, diff_method)
-
-    # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]²
-    xs, ws = _gausslegendre(FP, rule.n)
-    wws = Iterators.product(ws, ws)
-    xxs = Iterators.product(xs, xs)
-
-    # Normalize the Plane's orthogonal vectors
-    uu = Meshes.unormalize(plane.u)
-    uv = Meshes.unormalize(plane.v)
-    uplane = Meshes.Plane(plane.p, uu, uv)
-
-    # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞)
-    t(x) = x / (1 - x^2)
-    t′(x) = (1 + x^2) / (1 - x^2)^2
-
-    # Integrate f over the Plane
-    domainunits = _units(uplane(0, 0))
-    function integrand(((wi, wj), (xi, xj)))
-        wi * wj * f(uplane(t(xi), t(xj))) * t′(xi) * t′(xj) * domainunits^2
-    end
-    return sum(integrand, zip(wws, xxs))
+        rule::IntegrationRule;
+        kwargs...
+)
+    # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]²
+    param_plane = _ParametricGeometry(_parametric(plane), Meshes.paramdim(plane))
+
+    # Integrate the _ParametricGeometry using the standard methods
+    return _integral(f, param_plane, rule; kwargs...)
 end
 
-function integral(
-        f,
-        plane::Meshes.Plane,
-        rule::GaussKronrod;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Plane, diff_method)
-
-    # Normalize the Plane's orthogonal vectors
-    uu = Meshes.unormalize(plane.u)
-    uv = Meshes.unormalize(plane.v)
-    uplane = Meshes.Plane(plane.p, uu, uv)
-
-    # Integrate f over the Plane
-    domainunits = _units(uplane(0, 0))^2
-    integrand(u, v) = f(uplane(u, v)) * domainunits
-    inner∫(v) = QuadGK.quadgk(u -> integrand(u, v), FP(-Inf), FP(Inf); rule.kwargs...)[1]
-    return QuadGK.quadgk(inner∫, FP(-Inf), FP(Inf); rule.kwargs...)[1]
+############################################################################################
+#                                       Parametric
+############################################################################################
+
+# Map argument domain from [0, 1]² to (-∞, ∞)² for (::Plane)(t1, t2)
+function _parametric(plane::Meshes.Plane)
+    # [-1, 1] ↦ (-∞, ∞)
+    f1(t) = t / (1 - t^2)
+    # [0, 1] ↦ [-1, 1]
+    f2(t) = 2t - 1
+    # Compose the two transforms
+    f = f1 ∘ f2
+    return (t1, t2) -> plane(f(t1), f(t2))
 end
-
-function integral(
-        f,
-        plane::Meshes.Plane,
-        rule::HAdaptiveCubature;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Plane, diff_method)
-
-    # Normalize the Plane's orthogonal vectors
-    uu = Meshes.unormalize(plane.u)
-    uv = Meshes.unormalize(plane.v)
-    uplane = Meshes.Plane(plane.p, uu, uv)
-
-    # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞)
-    t(x) = x / (1 - x^2)
-    t′(x) = (1 + x^2) / (1 - x^2)^2
-
-    # Integrate f over the Plane
-    domainunits = _units(uplane(0, 0))
-    function integrand(xs)
-        f(uplane(t(xs[1]), t(xs[2]))) * t′(xs[1]) * t′(xs[2]) * domainunits^2
-    end
-
-    # HCubature doesn't support functions that output Unitful Quantity types
-    # Establish the units that are output by f
-    testpoint_parametriccoord = zeros(FP, 2)
-    integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
-    # Create a wrapper that returns only the value component in those units
-    uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
-    # Integrate only the unitless values
-    a = .-_ones(FP, 2)
-    b = _ones(FP, 2)
-    value = HCubature.hcubature(uintegrand, a, b; rule.kwargs...)[1]
-
-    # Reapply units
-    return value .* integrandunits
-end
-
-################################################################################
-#                               jacobian
-################################################################################
-
-_has_analytical(::Type{T}) where {T <: Meshes.Plane} = true

src/specializations/Ray.jl

@@ -11,85 +11,24 @@
 function integral(
         f,
         ray::Meshes.Ray,
-        rule::GaussLegendre;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Ray, diff_method)
-
-    # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
-    xs, ws = _gausslegendre(FP, rule.n)
-
-    # Normalize the Ray s.t. ray(t) is distance t from origin point
-    ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v))
-
-    # Domain transformation: x ∈ [-1,1] ↦ t ∈ [0,∞)
-    t₁(x) = (1 // 2) * x + (1 // 2)
-    t₁′(x) = (1 // 2)
-    t₂(x) = x / (1 - x^2)
-    t₂′(x) = (1 + x^2) / (1 - x^2)^2
-    t = t₂ ∘ t₁
-    t′(x) = t₂′(t₁(x)) * t₁′(x)
-
-    # Integrate f along the Ray
-    domainunits = _units(ray(0))
-    integrand(x) = f(ray(t(x))) * t′(x) * domainunits
-    return sum(w .* integrand(x) for (w, x) in zip(ws, xs))
-end
-
-function integral(
-        f,
-        ray::Meshes.Ray,
-        rule::GaussKronrod;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Ray, diff_method)
-
-    # Normalize the Ray s.t. ray(t) is distance t from origin point
-    ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v))
-
-    # Integrate f along the Ray
-    domainunits = _units(ray(0))
-    integrand(t) = f(ray(t)) * domainunits
-    return QuadGK.quadgk(integrand, zero(FP), FP(Inf); rule.kwargs...)[1]
-end
-
-function integral(
-        f,
-        ray::Meshes.Ray,
-        rule::HAdaptiveCubature;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Ray, diff_method)
-
-    # Normalize the Ray s.t. ray(t) is distance t from origin point
-    ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v))
-
-    # Domain transformation: x ∈ [0,1] ↦ t ∈ [0,∞)
-    t(x) = x / (1 - x^2)
-    t′(x) = (1 + x^2) / (1 - x^2)^2
-
-    # Integrate f along the Ray
-    domainunits = _units(ray(0))
-    integrand(xs) = f(ray(t(xs[1]))) * t′(xs[1]) * domainunits
-
-    # 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(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
-    # 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;
+        kwargs...
+)
+    # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]
+    param_ray = _ParametricGeometry(_parametric(ray), Meshes.paramdim(ray))
+
+    # Integrate the _ParametricGeometry using the standard methods
+    return _integral(f, param_ray, rule; kwargs...)
 end
 
 ################################################################################
-#                               jacobian
+#                                 Parametric
 ################################################################################
 
-_has_analytical(::Type{T}) where {T <: Meshes.Ray} = true
+# Map argument domain from [0, 1] to [0, ∞) for (::Ray)(t)
+# f(t) = t / 1 - t^2)
+# f'(t) = (t^2 + 1) / (1 - t^2)^2
+function _parametric(ray::Meshes.Ray)
+    f(t) = t / (1 - t^2)
+    return t -> ray(f(t))
+end