From 93b1971db5d4c9f98ea7b0676d4ee039674e5b2e Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 14:53:29 -0500
Subject: [PATCH 01/33] Implement _ParametricGeometry for Line

---
 src/specializations/Line.jl | 89 +++++++------------------------------
 1 file changed, 16 insertions(+), 73 deletions(-)

diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl
index 0c7da2c1..ab88b86b 100644
--- a/src/specializations/Line.jl
+++ b/src/specializations/Line.jl
@@ -10,81 +10,24 @@
 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...
+)
+    paramfunction(t) = _parametric(line, t)
+    param_line = _ParametricGeometry(paramfunction, 1)
+    return _integral(f, param_line, rule; kwargs...)
 end
 
 ################################################################################
-#                               jacobian
+#                              Parametric
 ################################################################################
 
-_has_analytical(::Type{T}) where {T <: Meshes.Line} = true
+# Map [0, 1] ↦ (-∞, ∞)
+function _parametric(line::Meshes.Line, t)
+    # [-1, 1] ↦ (-∞, ∞)
+    f1(t) = t / (1 - t^2)
+    # [0, 1] ↦ [-1, 1]
+    f2(t) = 2t - 1
+    # Compose the two transforms
+    return line((f1 ∘ f2)(t))
+end

From 8f3502b430fb766105c48aaff199014c79a1818a Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 16:18:12 -0500
Subject: [PATCH 02/33] Implement _ParametricGeometry for Plane

---
 src/specializations/Plane.jl | 104 +++++------------------------------
 1 file changed, 15 insertions(+), 89 deletions(-)

diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl
index ee2706cb..77c30584 100644
--- a/src/specializations/Plane.jl
+++ b/src/specializations/Plane.jl
@@ -11,96 +11,22 @@
 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))
-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]
+        rule::IntegrationRule;
+        kwargs...
+)
+    paramfunction(t1, t2) = _parametric(plane, t1, t2)
+    param_plane = _ParametricGeometry(paramfunction, 2)
+    return _integral(f, param_plane, rule; kwargs...)
 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
+############################################################################################
+#                                       Parametric
+############################################################################################
 
-    # 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
+# Map [0, 1]² ↦ (-∞, ∞)²
+function _parametric(plane::Meshes.Plane, t1, t2)
+    f1(t) = t / (1 - t^2)
+    f2(t) = 2t - 1
+    f = f1 ∘ f2
+    return plane(f(t1), f(t2))
 end
-
-################################################################################
-#                               jacobian
-################################################################################
-
-_has_analytical(::Type{T}) where {T <: Meshes.Plane} = true

From 201ae62787156126d98f8cc3d0a59ada0dc58030 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 16:18:18 -0500
Subject: [PATCH 03/33] Implement _ParametricGeometry for Ray

---
 src/specializations/Ray.jl | 91 ++++++--------------------------------
 1 file changed, 13 insertions(+), 78 deletions(-)

diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index eede812e..8afc9fc2 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -11,85 +11,20 @@
 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]
+        rule::IntegrationRule;
+        kwargs...
+)
+    paramfunction(t) = _parametric(ray, t)
+    param_ray = _ParametricGeometry(paramfunction, 1)
+    return _integral(f, param_ray, rule; kwargs...)
 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
+############################################################################################
+#                                     Parametric
+############################################################################################
 
-    # 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
+# Map [0, 1] ↦ [0, ∞)
+function _parametric(ray::Meshes.Ray, t)
+    f(t) = t / (1 - t^2)
+    return ray(f(t))
 end
-
-################################################################################
-#                               jacobian
-################################################################################
-
-_has_analytical(::Type{T}) where {T <: Meshes.Ray} = true

From cdac6ec449994f3ac07fdae5e058b9c0aab02135 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 16:21:59 -0500
Subject: [PATCH 04/33] Update and reorganize analytical tests

---
 test/utils.jl | 29 +++++++++++++----------------
 1 file changed, 13 insertions(+), 16 deletions(-)

diff --git a/test/utils.jl b/test/utils.jl
index 69a8fde6..f02dd455 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -22,30 +22,27 @@ end
 @testitem "DifferentiationMethod" setup=[Setup] begin
     using MeshIntegrals: _has_analytical, _default_method, _guarantee_analytical
 
-    # _has_analytical of instances
-    triangle = Triangle(Point(0, 0, 0), Point(0, 1, 0), Point(1, 0, 0))
-    @test _has_analytical(triangle) == true
+    # Test geometries
     sphere = Sphere(Point(0, 0, 0), 1.0)
-    @test _has_analytical(sphere) == false
+    triangle = Triangle(Point(0, 0, 0), Point(0, 1, 0), Point(1, 0, 0))
 
-    # _default_method
-    @test _default_method(Meshes.Triangle) isa Analytical
-    @test _default_method(triangle) isa Analytical
-    @test _default_method(Meshes.Sphere) isa FiniteDifference
-    @test _default_method(sphere) isa FiniteDifference
+    # _has_analytical of instances
+    @test _has_analytical(sphere) == false
+    @test _has_analytical(triangle) == true
 
     # _has_analytical of types
-    @test _has_analytical(Meshes.BezierCurve) == false
-    @test _has_analytical(Meshes.Line) == true
-    @test _has_analytical(Meshes.Plane) == true
-    @test _has_analytical(Meshes.Ray) == true
     @test _has_analytical(Meshes.Sphere) == false
-    @test _has_analytical(Meshes.Tetrahedron) == false
     @test _has_analytical(Meshes.Triangle) == true
 
+    # _default_method
+    @test _default_method(Meshes.Sphere) isa FiniteDifference
+    @test _default_method(sphere) isa FiniteDifference
+    @test _default_method(Meshes.Triangle) isa Analytical
+    @test _default_method(triangle) isa Analytical
+
     # _guarantee_analytical
-    @test _guarantee_analytical(Meshes.Line, Analytical()) === nothing
-    @test_throws "Analytical" _guarantee_analytical(Meshes.Line, FiniteDifference())
+    @test _guarantee_analytical(Meshes.Triangle, Analytical()) === nothing
+    @test_throws "Analytical" _guarantee_analytical(Meshes.Triangle, FiniteDifference())
 
     # FiniteDifference
     @test FiniteDifference().ε ≈ 1e-6

From f4518c89713ac5d717d08c7f20a8a3a35021d2ff Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 16:40:27 -0500
Subject: [PATCH 05/33] Add comments

---
 src/specializations/BezierCurve.jl | 1 +
 src/specializations/Line.jl        | 6 ++++--
 src/specializations/Plane.jl       | 9 +++++++--
 src/specializations/Ray.jl         | 6 ++++--
 4 files changed, 16 insertions(+), 6 deletions(-)

diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index 01b0f24d..d062c403 100644
--- a/src/specializations/BezierCurve.jl
+++ b/src/specializations/BezierCurve.jl
@@ -50,6 +50,7 @@ end
 #                              Parametric
 ################################################################################
 
+# Wrap (::BezierCurve)(t, ::BezierEvalMethod) into f(t) by embedding second argument
 function _parametric(curve::Meshes.BezierCurve, t, alg::Meshes.BezierEvalMethod)
     return curve(t, alg)
 end
diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl
index ab88b86b..b1473b49 100644
--- a/src/specializations/Line.jl
+++ b/src/specializations/Line.jl
@@ -13,8 +13,10 @@ function integral(
         rule::IntegrationRule;
         kwargs...
 )
+    # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]
     paramfunction(t) = _parametric(line, t)
-    param_line = _ParametricGeometry(paramfunction, 1)
+
+    # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_line, rule; kwargs...)
 end
 
@@ -22,7 +24,7 @@ end
 #                              Parametric
 ################################################################################
 
-# Map [0, 1] ↦ (-∞, ∞)
+# Map argument domain from [0, 1] to (-∞, ∞) for (::Line)(t)
 function _parametric(line::Meshes.Line, t)
     # [-1, 1] ↦ (-∞, ∞)
     f1(t) = t / (1 - t^2)
diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl
index 77c30584..5825aa45 100644
--- a/src/specializations/Plane.jl
+++ b/src/specializations/Plane.jl
@@ -14,8 +14,10 @@ function integral(
         rule::IntegrationRule;
         kwargs...
 )
+    # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]²
     paramfunction(t1, t2) = _parametric(plane, t1, t2)
-    param_plane = _ParametricGeometry(paramfunction, 2)
+
+    # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_plane, rule; kwargs...)
 end
 
@@ -23,10 +25,13 @@ end
 #                                       Parametric
 ############################################################################################
 
-# Map [0, 1]² ↦ (-∞, ∞)²
+# Map argument domain from [0, 1]² to (-∞, ∞)² for (::Plane)(t1, t2)
 function _parametric(plane::Meshes.Plane, t1, t2)
+    # [-1, 1] ↦ (-∞, ∞)
     f1(t) = t / (1 - t^2)
+    # [0, 1] ↦ [-1, 1]
     f2(t) = 2t - 1
+    # Compose the two transforms
     f = f1 ∘ f2
     return plane(f(t1), f(t2))
 end
diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index 8afc9fc2..736e3a05 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -14,8 +14,10 @@ function integral(
         rule::IntegrationRule;
         kwargs...
 )
+    # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]
     paramfunction(t) = _parametric(ray, t)
-    param_ray = _ParametricGeometry(paramfunction, 1)
+
+    # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_ray, rule; kwargs...)
 end
 
@@ -23,7 +25,7 @@ end
 #                                     Parametric
 ############################################################################################
 
-# Map [0, 1] ↦ [0, ∞)
+# Map argument domain from [0, 1] to [0, ∞) for (::Ray)(t)
 function _parametric(ray::Meshes.Ray, t)
     f(t) = t / (1 - t^2)
     return ray(f(t))

From 3ce948c71b0cca08fd8331c1241d78d0c607a636 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 16:40:49 -0500
Subject: [PATCH 06/33] Use paramdim for clarity

---
 src/specializations/Line.jl  | 1 +
 src/specializations/Plane.jl | 1 +
 src/specializations/Ray.jl   | 1 +
 3 files changed, 3 insertions(+)

diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl
index b1473b49..9dcfae4e 100644
--- a/src/specializations/Line.jl
+++ b/src/specializations/Line.jl
@@ -15,6 +15,7 @@ function integral(
 )
     # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]
     paramfunction(t) = _parametric(line, t)
+    param_line = _ParametricGeometry(paramfunction, Meshes.paramdim(line))
 
     # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_line, rule; kwargs...)
diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl
index 5825aa45..22673957 100644
--- a/src/specializations/Plane.jl
+++ b/src/specializations/Plane.jl
@@ -16,6 +16,7 @@ function integral(
 )
     # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]²
     paramfunction(t1, t2) = _parametric(plane, t1, t2)
+    param_plane = _ParametricGeometry(paramfunction, Meshes.paramdim(plane))
 
     # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_plane, rule; kwargs...)
diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index 736e3a05..0446fa93 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -16,6 +16,7 @@ function integral(
 )
     # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]
     paramfunction(t) = _parametric(ray, t)
+    param_ray = _ParametricGeometry(paramfunction, Meshes.paramdim(ray))
 
     # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_ray, rule; kwargs...)

From 7ce264e04471ab5355c43b80ec26a1d6a746d317 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 16:45:38 -0500
Subject: [PATCH 07/33] Test _parametric -> anonymous function

---
 src/specializations/Ray.jl | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index 0446fa93..4fbd0e83 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -15,8 +15,7 @@ function integral(
         kwargs...
 )
     # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]
-    paramfunction(t) = _parametric(ray, t)
-    param_ray = _ParametricGeometry(paramfunction, Meshes.paramdim(ray))
+    param_ray = _ParametricGeometry(_parametric(ray), Meshes.paramdim(ray))
 
     # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_ray, rule; kwargs...)
@@ -27,7 +26,7 @@ end
 ############################################################################################
 
 # Map argument domain from [0, 1] to [0, ∞) for (::Ray)(t)
-function _parametric(ray::Meshes.Ray, t)
+function _parametric(ray::Meshes.Ray)
     f(t) = t / (1 - t^2)
-    return ray(f(t))
+    return t -> ray(f(t))
 end

From 8e9c6c87b6a9c6e945f40da7f2e01d1d34b10c34 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 18:08:51 -0500
Subject: [PATCH 08/33] Update _parametric API to return an anonymous function

---
 src/specializations/BezierCurve.jl         | 7 +++----
 src/specializations/Line.jl                | 7 +++----
 src/specializations/Plane.jl               | 7 +++----
 src/specializations/_ParametricGeometry.jl | 2 +-
 4 files changed, 10 insertions(+), 13 deletions(-)

diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index d062c403..a79f5dd3 100644
--- a/src/specializations/BezierCurve.jl
+++ b/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...)
@@ -51,6 +50,6 @@ end
 ################################################################################
 
 # Wrap (::BezierCurve)(t, ::BezierEvalMethod) into f(t) by embedding second argument
-function _parametric(curve::Meshes.BezierCurve, t, alg::Meshes.BezierEvalMethod)
-    return curve(t, alg)
+function _parametric(curve::Meshes.BezierCurve, alg::Meshes.BezierEvalMethod)
+    return t -> curve(t, alg)
 end
diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl
index 9dcfae4e..e099928c 100644
--- a/src/specializations/Line.jl
+++ b/src/specializations/Line.jl
@@ -14,8 +14,7 @@ function integral(
         kwargs...
 )
     # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]
-    paramfunction(t) = _parametric(line, t)
-    param_line = _ParametricGeometry(paramfunction, Meshes.paramdim(line))
+    param_line = _ParametricGeometry(_parametric(line), Meshes.paramdim(line))
 
     # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_line, rule; kwargs...)
@@ -26,11 +25,11 @@ end
 ################################################################################
 
 # Map argument domain from [0, 1] to (-∞, ∞) for (::Line)(t)
-function _parametric(line::Meshes.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 line((f1 ∘ f2)(t))
+    return t -> line((f1 ∘ f2)(t))
 end
diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl
index 22673957..df7c4e2d 100644
--- a/src/specializations/Plane.jl
+++ b/src/specializations/Plane.jl
@@ -15,8 +15,7 @@ function integral(
         kwargs...
 )
     # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]²
-    paramfunction(t1, t2) = _parametric(plane, t1, t2)
-    param_plane = _ParametricGeometry(paramfunction, Meshes.paramdim(plane))
+    param_plane = _ParametricGeometry(_parametric(plane), Meshes.paramdim(plane))
 
     # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_plane, rule; kwargs...)
@@ -27,12 +26,12 @@ end
 ############################################################################################
 
 # Map argument domain from [0, 1]² to (-∞, ∞)² for (::Plane)(t1, t2)
-function _parametric(plane::Meshes.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 plane(f(t1), f(t2))
+    return (t1, t2) -> plane(f(t1), f(t2))
 end
diff --git a/src/specializations/_ParametricGeometry.jl b/src/specializations/_ParametricGeometry.jl
index 16904d27..bdf0f6c9 100644
--- a/src/specializations/_ParametricGeometry.jl
+++ b/src/specializations/_ParametricGeometry.jl
@@ -55,7 +55,7 @@ end
 Meshes.paramdim(::_ParametricGeometry{M, C, F, Dim}) where {M, C, F, Dim} = Dim
 
 """
-    _parametric(geometry::G, ts...) where {G <: Meshes.Geometry}
+    _parametric(geometry::G) where {G <: Meshes.Geometry} -> Function
 
 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.

From 77f39f460274c0575ae1805618f6d963a1e133a8 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 18:22:19 -0500
Subject: [PATCH 09/33] Remove obsolete N, enable Tetrahedron test

---
 benchmark/benchmarks.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index f8f70ee9..4b3c44bf 100644
--- a/benchmark/benchmarks.jl
+++ b/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)
     s
 end
 

From 0b3739adc592ef00a2c961fb1c7cb265bec18d13 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 18:26:26 -0500
Subject: [PATCH 10/33] Remove hanging comma

---
 benchmark/benchmarks.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 4b3c44bf..eb10808f 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -26,7 +26,7 @@ geometries = (
 
 SUITE["Integrals"] = let s = BenchmarkGroup()
     for (int, rule, geometry) in Iterators.product(integrands, rules, geometries)
-        n1 = geometry.name,
+        n1 = geometry.name
         n2 = "$(int.name) $(rule.name)"
         s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule)
     end

From 2e06577b2c29b9216e410d480aa482cbd341d2a4 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 18:31:58 -0500
Subject: [PATCH 11/33] Update field name

---
 benchmark/benchmarks.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index eb10808f..b528be0c 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -65,7 +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)
+    s["Tetrahedron"] = @benchmarkable integral($spec.f, $spec.g.tetrahedron, $spec.rule_gl)
     s
 end
 

From 5f404f3ccec0beae92f98af7aa3827eabc05183b Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 18:46:04 -0500
Subject: [PATCH 12/33] Add analytical derivative in comment

---
 src/specializations/Ray.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index 4fbd0e83..825f2a71 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -21,11 +21,13 @@ function integral(
     return _integral(f, param_ray, rule; kwargs...)
 end
 
-############################################################################################
-#                                     Parametric
-############################################################################################
+################################################################################
+#                                 Parametric
+################################################################################
 
 # 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))

From dce7c3a5146bb337eb7677eb485639119700a359 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 19:12:26 -0500
Subject: [PATCH 13/33] Update _parametric for Triangle and Tetrahedron

---
 src/specializations/Tetrahedron.jl | 30 ++++++++++++++++--------------
 src/specializations/Triangle.jl    | 25 ++++++++++++++-----------
 2 files changed, 30 insertions(+), 25 deletions(-)

diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl
index 2f6fdece..d3572c53 100644
--- a/src/specializations/Tetrahedron.jl
+++ b/src/specializations/Tetrahedron.jl
@@ -15,9 +15,8 @@ function integral(
         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))
+    # Generate a _ParametricGeometry whose parametric function domain spans [0,1]³
+    tetra = _ParametricGeometry(_parametric(tetrahedron), Meshes.paramdim(tetrahedron))
 
     # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, tetra, rule; kwargs...)
@@ -27,17 +26,20 @@ end
 #                               Parametric
 ################################################################################
 
-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
+function _parametric(tetrahedron::Meshes.Tetrahedron)
+    function f(t1, t2, t3)
+        if any(Iterators.map(n -> (n < 0) || (n > 1), (t1, t2, t3)))
+            msg = "tetrahedron(t1, t2, t3) is not defined for t outside [0, 1]."
+            throw(DomainError((t1, t2, 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)
+        # Take a triangular cross-section at t3
+        a = tetrahedron(t3, 0, 0)
+        b = tetrahedron(0, t3, 0)
+        c = tetrahedron(0, 0, t3)
+        cross_section = _parametric(Meshes.Triangle(a, b, c))
 
-    return _parametric(cross_section, t1, t2)
+        return cross_section(t1, t2)
+    end
+    return f
 end
diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index 2c70f400..69f4cb7b 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -111,16 +111,19 @@ _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))
+function _parametric(triangle::Meshes.Triangle)
+    function f(t1, t2)
+        if any(Iterators.map(n -> (n < 0) || (n > 1), (t1, t2)))
+            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
-
-    # Form a line segment between sides
-    a = triangle(0, t2)
-    b = triangle(t2, 0)
-    segment = Meshes.Segment(a, b)
-
-    return segment(t1)
+    return f
 end

From 63e85fefbcf3d058a3b95ff011ef7f2a5411e642 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Sat, 30 Nov 2024 19:29:19 -0500
Subject: [PATCH 14/33] Update tests

---
 test/utils.jl | 18 +++++++-----------
 1 file changed, 7 insertions(+), 11 deletions(-)

diff --git a/test/utils.jl b/test/utils.jl
index f02dd455..afd8d9f0 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -55,18 +55,14 @@ end
     pt_w = Point(-7, 0, 0)
     pt_e = Point(8, 0, 0)
     ẑ = Vec(0, 0, 1)
-    triangle = Triangle(pt_n, pt_w, pt_e)
-    tetrahedron = Tetrahedron(pt_n, pt_w, pt_e, pt_n + ẑ)
+    triangle = _parametric(Triangle(pt_n, pt_w, pt_e))
+    tetrahedron = _parametric(Tetrahedron(pt_n, pt_w, pt_e, pt_n + ẑ))
 
-    # Ensure error is thrown for an out-of-bounds coordinate
-    for ts in [(-0.5, 0.5), (0.5, -0.5), (1.5, 0.5), (0.5, 1.5)]
-        @test_throws "not defined" _parametric(triangle, ts...)
-        # Tetrahedron forwards t1 and t2 to _parametric(::Triangle, ts...)
-        @test_throws "not defined" _parametric(tetrahedron, ts..., 0)
+    # Ensure error is thrown for out-of-bounds coordinate
+    for ts in [(-1, 0), (0, -1), (2, 0), (0, 2)]
+        @test_throws "not defined" triangle(ts...)
     end
-
-    # Ensue error is thrown for an out-of-bounds third coordinate
-    for t3 in [-0.5, 1.5]
-        @test_throws "not defined" _parametric(tetrahedron, 0, 0, t3)
+    for ts in [(-1, 0, 0), (0, -1, 0), (0, 0, -1), (2, 0, 0), (0, 2, 0), (0, 0, 2)]
+        @test_throws "not defined" tetrahedron(ts...)
     end
 end

From b5c6886e0dcae94cf172a9832f0ef706f47cde6e Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Dec 2024 19:27:27 +0100
Subject: [PATCH 15/33] put parametrized to .typos.toml

---
 .typos.toml | 2 ++
 1 file changed, 2 insertions(+)
 create mode 100644 .typos.toml

diff --git a/.typos.toml b/.typos.toml
new file mode 100644
index 00000000..0f6b1ff9
--- /dev/null
+++ b/.typos.toml
@@ -0,0 +1,2 @@
+[default.extend-words]
+parametrized = "parametrized"

From 4217d11cf008c89f4b8f308215972e5398e8c12c Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 20:10:16 -0500
Subject: [PATCH 16/33] Implement _ParametricGeometry for Triangle

---
 src/specializations/Triangle.jl | 107 ++++----------------------------
 1 file changed, 11 insertions(+), 96 deletions(-)

diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index 69f4cb7b..c4c55416 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -10,103 +10,18 @@
 ################################################################################
 
 function integral(
-        f,
-        triangle::Meshes.Triangle,
-        rule::GaussLegendre;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Triangle, diff_method)
-
-    # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2
-    xs, ws = _gausslegendre(FP, rule.n)
-    wws = Iterators.product(ws, ws)
-    xxs = Iterators.product(xs, xs)
-
-    # Domain transformations:
-    #   xᵢ [-1,1] ↦ R [0,1]
-    #   xⱼ [-1,1] ↦ φ [0,π/2]
-    uR(xᵢ) = (1 // 2) * (xᵢ + 1)
-    uφ(xⱼ) = (T(π) / 4) * (xⱼ + 1)
-
-    # Integrate the Barycentric triangle by transforming it into polar coordinates
-    #   with a modified radius
-    #     R = r ( sinφ + cosφ )
-    #   s.t. integration bounds become rectangular
-    #     R ∈ [0, 1] and φ ∈ [0, π/2]
-    function integrand(((wᵢ, wⱼ), (xᵢ, xⱼ)))
-        R = uR(xᵢ)
-        φ = uφ(xⱼ)
-        a, b = sincos(φ)
-        u = R * (1 - a / (a + b))
-        v = R * (1 - b / (a + b))
-        return wᵢ * wⱼ * f(triangle(u, v)) * R / (a + b)^2
-    end
-
-    # Calculate 2D Gauss-Legendre integral over modified-polar-Barycentric coordinates
-    # Apply a linear domain-correction factor
-    return FP(π / 4) * Meshes.area(triangle) .* sum(integrand, zip(wws, xxs))
+    f,
+    triangle::Meshes.Triangle,
+    rule::IntegrationRule;
+    kwargs...
+)
+    # Generate a _ParametricGeometry whose parametric function domain spans [0,1]²
+    param_triangle = _ParametricGeometry(_parametric(triangle), 2)
+
+    # Integrate the _ParametricGeometry using the standard methods
+    return _integral(f, param_triangle, rule; kwargs...)
 end
 
-function integral(
-        f,
-        triangle::Meshes.Triangle,
-        rule::GaussKronrod;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Triangle, diff_method)
-
-    # Integrate the Barycentric triangle in (u,v)-space: (0,0), (0,1), (1,0)
-    #   i.e. \int_{0}^{1} \int_{0}^{1-u} f(u,v) dv du
-    ∫u(u) = QuadGK.quadgk(v -> f(triangle(u, v)), zero(FP), FP(1 - u); rule.kwargs...)[1]
-    ∫ = QuadGK.quadgk(∫u, zero(FP), one(FP); rule.kwargs...)[1]
-
-    # Apply a linear domain-correction factor 0.5 ↦ area(triangle)
-    return 2 * Meshes.area(triangle) .* ∫
-end
-
-function integral(
-        f,
-        triangle::Meshes.Triangle,
-        rule::HAdaptiveCubature;
-        diff_method::DM = Analytical(),
-        FP::Type{T} = Float64
-) where {DM <: DifferentiationMethod, T <: AbstractFloat}
-    _guarantee_analytical(Meshes.Triangle, diff_method)
-
-    # Integrate the Barycentric triangle by transforming it into polar coordinates
-    #   with a modified radius
-    #     R = r ( sinφ + cosφ )
-    #   s.t. integration bounds become rectangular
-    #     R ∈ [0, 1] and φ ∈ [0, π/2]
-    function integrand(Rφ)
-        R, φ = Rφ
-        a, b = sincos(φ)
-        u = R * (1 - a / (a + b))
-        v = R * (1 - b / (a + b))
-        return f(triangle(u, v)) * R / (a + b)^2
-    end
-
-    # 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) * integrandunits .* ∫
-end
-
-################################################################################
-#                               jacobian
-################################################################################
-
-_has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true
-
 ################################################################################
 #                              Parametric
 ################################################################################
@@ -114,7 +29,7 @@ _has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true
 function _parametric(triangle::Meshes.Triangle)
     function f(t1, t2)
         if any(Iterators.map(n -> (n < 0) || (n > 1), (t1, t2)))
-            msg = "triangle(t1, t2) is not defined for (t1, t2) outside [0, 1]^2."
+            msg = "triangle(t1, t2) is not defined for (t1, t2) outside [0, 1]²."
             throw(DomainError((t1, t2), msg))
         end
 

From 8b8ff2884992eec1cf97d6d509c3dc4991016609 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 20:10:38 -0500
Subject: [PATCH 17/33] Temp disable some Analytical tests

---
 test/utils.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/utils.jl b/test/utils.jl
index afd8d9f0..323767c9 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -22,6 +22,7 @@ end
 @testitem "DifferentiationMethod" setup=[Setup] begin
     using MeshIntegrals: _has_analytical, _default_method, _guarantee_analytical
 
+    #=
     # Test geometries
     sphere = Sphere(Point(0, 0, 0), 1.0)
     triangle = Triangle(Point(0, 0, 0), Point(0, 1, 0), Point(1, 0, 0))
@@ -33,6 +34,7 @@ end
     # _has_analytical of types
     @test _has_analytical(Meshes.Sphere) == false
     @test _has_analytical(Meshes.Triangle) == true
+    =#
 
     # _default_method
     @test _default_method(Meshes.Sphere) isa FiniteDifference
@@ -40,9 +42,11 @@ end
     @test _default_method(Meshes.Triangle) isa Analytical
     @test _default_method(triangle) isa Analytical
 
+    #=
     # _guarantee_analytical
     @test _guarantee_analytical(Meshes.Triangle, Analytical()) === nothing
     @test_throws "Analytical" _guarantee_analytical(Meshes.Triangle, FiniteDifference())
+    =#
 
     # FiniteDifference
     @test FiniteDifference().ε ≈ 1e-6

From 1cac8308c6f039b66ed9bca4358e034932f87646 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 20:44:25 -0500
Subject: [PATCH 18/33] Temp disable more tests

---
 test/utils.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/utils.jl b/test/utils.jl
index 323767c9..8fbaad0f 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -39,8 +39,8 @@ end
     # _default_method
     @test _default_method(Meshes.Sphere) isa FiniteDifference
     @test _default_method(sphere) isa FiniteDifference
-    @test _default_method(Meshes.Triangle) isa Analytical
-    @test _default_method(triangle) isa Analytical
+    #@test _default_method(Meshes.Triangle) isa Analytical
+    #@test _default_method(triangle) isa Analytical
 
     #=
     # _guarantee_analytical

From 044a20de961897cb94e1689a815464bacc3fa1d9 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 20:45:32 -0500
Subject: [PATCH 19/33] Test disabling bound check

---
 src/specializations/Triangle.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index c4c55416..bbd1e4f6 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -28,10 +28,12 @@ end
 
 function _parametric(triangle::Meshes.Triangle)
     function f(t1, t2)
+        #=
         if any(Iterators.map(n -> (n < 0) || (n > 1), (t1, t2)))
             msg = "triangle(t1, t2) is not defined for (t1, t2) outside [0, 1]²."
             throw(DomainError((t1, t2), msg))
         end
+        =#
 
         # Form a line segment between sides
         a = triangle(0, t2)

From 15710a64ad03d890a36c1c6506af0361dfad14ec Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 21:01:34 -0500
Subject: [PATCH 20/33] Restore bounds check, try new formulation

---
 src/specializations/Triangle.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index bbd1e4f6..a9af8e7a 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -28,19 +28,20 @@ end
 
 function _parametric(triangle::Meshes.Triangle)
     function f(t1, t2)
-        #=
         if any(Iterators.map(n -> (n < 0) || (n > 1), (t1, t2)))
             msg = "triangle(t1, t2) is not defined for (t1, t2) outside [0, 1]²."
             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)
+        =#
+        return triangle(t1 * t2, t2 - (t1 * t2))
     end
     return f
 end

From 7d3c83e437935fd3e32ed096f665f21e6be16b57 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 21:16:36 -0500
Subject: [PATCH 21/33] Remove old commented-out code

---
 src/specializations/Triangle.jl | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index a9af8e7a..0d0fc133 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -33,14 +33,6 @@ function _parametric(triangle::Meshes.Triangle)
             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)
-        =#
         return triangle(t1 * t2, t2 - (t1 * t2))
     end
     return f

From dbdffedb26b9986cc30172470bd7bc4b9efce9fb Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 21:20:04 -0500
Subject: [PATCH 22/33] Add comments

---
 src/specializations/Tetrahedron.jl | 1 +
 src/specializations/Triangle.jl    | 1 +
 2 files changed, 2 insertions(+)

diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl
index d3572c53..1aedfaf5 100644
--- a/src/specializations/Tetrahedron.jl
+++ b/src/specializations/Tetrahedron.jl
@@ -26,6 +26,7 @@ end
 #                               Parametric
 ################################################################################
 
+# Map argument domain from [0, 1]³ to Barycentric domain for (::Tetrahedron)(t1, t2, t3)
 function _parametric(tetrahedron::Meshes.Tetrahedron)
     function f(t1, t2, t3)
         if any(Iterators.map(n -> (n < 0) || (n > 1), (t1, t2, t3)))
diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index 0d0fc133..f8805405 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -26,6 +26,7 @@ end
 #                              Parametric
 ################################################################################
 
+# Map argument domain from [0, 1]² to Barycentric domain for (::Triangle)(t1, t2)
 function _parametric(triangle::Meshes.Triangle)
     function f(t1, t2)
         if any(Iterators.map(n -> (n < 0) || (n > 1), (t1, t2)))

From bc4e591e319693a61dc41e6e2117c3b575097310 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 21:20:25 -0500
Subject: [PATCH 23/33] Remove obsolete tests

---
 test/utils.jl | 18 ------------------
 1 file changed, 18 deletions(-)

diff --git a/test/utils.jl b/test/utils.jl
index 8fbaad0f..ad336332 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -22,31 +22,13 @@ end
 @testitem "DifferentiationMethod" setup=[Setup] begin
     using MeshIntegrals: _has_analytical, _default_method, _guarantee_analytical
 
-    #=
     # Test geometries
     sphere = Sphere(Point(0, 0, 0), 1.0)
     triangle = Triangle(Point(0, 0, 0), Point(0, 1, 0), Point(1, 0, 0))
 
-    # _has_analytical of instances
-    @test _has_analytical(sphere) == false
-    @test _has_analytical(triangle) == true
-
-    # _has_analytical of types
-    @test _has_analytical(Meshes.Sphere) == false
-    @test _has_analytical(Meshes.Triangle) == true
-    =#
-
     # _default_method
     @test _default_method(Meshes.Sphere) isa FiniteDifference
     @test _default_method(sphere) isa FiniteDifference
-    #@test _default_method(Meshes.Triangle) isa Analytical
-    #@test _default_method(triangle) isa Analytical
-
-    #=
-    # _guarantee_analytical
-    @test _guarantee_analytical(Meshes.Triangle, Analytical()) === nothing
-    @test_throws "Analytical" _guarantee_analytical(Meshes.Triangle, FiniteDifference())
-    =#
 
     # FiniteDifference
     @test FiniteDifference().ε ≈ 1e-6

From 3178e8fb8f96bbc9584841b84e6c7abf02a1cb9a Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 21:21:43 -0500
Subject: [PATCH 24/33] Remove obsolete utils and update default diff_method

---
 src/utils.jl | 16 +---------------
 1 file changed, 1 insertion(+), 15 deletions(-)

diff --git a/src/utils.jl b/src/utils.jl
index e69a38ef..1e5292f8 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -26,25 +26,11 @@ _ones(N::Int) = _ones(Float64, N)
 #                           DifferentiationMethod
 ################################################################################
 
-# Throw an ArgumentError if Analytical() jacobian not defined for this type
-function _guarantee_analytical(
-        ::Type{G},
-        ::DifferentiationMethod
-) where {G <: Geometry}
-    throw(ArgumentError("$G geometries require kwarg diff_method = Analytical()"))
-end
-
-_guarantee_analytical(::Type{G}, ::Analytical) where {G <: Geometry} = nothing
-
-# Return whether a geometry type has jacobian methods defined
-_has_analytical(::Type{G}) where {G <: Geometry} = false
-_has_analytical(g::G) where {G <: Geometry} = _has_analytical(G)
-
 # Return the default DifferentiationMethod instance for a particular geometry type
 function _default_method(
         g::Type{G}
 ) where {G <: Geometry}
-    return _has_analytical(G) ? Analytical() : FiniteDifference()
+    return FiniteDifference()
 end
 
 # Return the default DifferentiationMethod instance for a particular geometry instance

From 96169010e766bc3824479d27e997ab225a32a631 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 21:34:58 -0500
Subject: [PATCH 25/33] Remove obsolete imports

---
 test/utils.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/utils.jl b/test/utils.jl
index ad336332..545b2540 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -20,7 +20,7 @@
 end
 
 @testitem "DifferentiationMethod" setup=[Setup] begin
-    using MeshIntegrals: _has_analytical, _default_method, _guarantee_analytical
+    using MeshIntegrals: _default_method
 
     # Test geometries
     sphere = Sphere(Point(0, 0, 0), 1.0)

From c21b433b15c6f7ecf19ebaae8a4848a475930690 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 21:35:39 -0500
Subject: [PATCH 26/33] Remove obsolete Analytical

---
 src/differentiation.jl | 16 ----------------
 1 file changed, 16 deletions(-)

diff --git a/src/differentiation.jl b/src/differentiation.jl
index 7ea44db0..5a5a0a6a 100644
--- a/src/differentiation.jl
+++ b/src/differentiation.jl
@@ -27,22 +27,6 @@ end
 FiniteDifference{T}() where {T <: AbstractFloat} = FiniteDifference{T}(T(1e-6))
 FiniteDifference() = FiniteDifference{Float64}()
 
-"""
-    Analytical()
-
-Use to specify use of analytically-derived solutions for calculating derivatives.
-These solutions are currently defined only for a subset of geometry types.
-
-# Supported Geometries:
-- `BezierCurve`
-- `Line`
-- `Plane`
-- `Ray`
-- `Tetrahedron`
-- `Triangle`
-"""
-struct Analytical <: DifferentiationMethod end
-
 # Future Support:
 #   struct AutoEnzyme <: DifferentiationMethod end
 #   struct AutoZygote <: DifferentiationMethod end

From 8a8f892c66539dcb949aa5e47cfa175034e97480 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 22:13:04 -0500
Subject: [PATCH 27/33] Remove obsolete export

---
 src/MeshIntegrals.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl
index 7b270c7b..cbc7260f 100644
--- a/src/MeshIntegrals.jl
+++ b/src/MeshIntegrals.jl
@@ -10,7 +10,7 @@ import QuadGK
 import Unitful
 
 include("differentiation.jl")
-export DifferentiationMethod, Analytical, FiniteDifference, jacobian
+export DifferentiationMethod, FiniteDifference, jacobian
 
 include("utils.jl")
 

From 6743a044a06b37f446ad25e202e82adabae5dde5 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 22:14:18 -0500
Subject: [PATCH 28/33] Formatting

---
 src/specializations/Ray.jl      | 2 --
 src/specializations/Triangle.jl | 8 ++++----
 2 files changed, 4 insertions(+), 6 deletions(-)

diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index 825f2a71..8c6e9e98 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -26,8 +26,6 @@ end
 ################################################################################
 
 # 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))
diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index f8805405..04308ed8 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -10,10 +10,10 @@
 ################################################################################
 
 function integral(
-    f,
-    triangle::Meshes.Triangle,
-    rule::IntegrationRule;
-    kwargs...
+        f,
+        triangle::Meshes.Triangle,
+        rule::IntegrationRule;
+        kwargs...
 )
     # Generate a _ParametricGeometry whose parametric function domain spans [0,1]²
     param_triangle = _ParametricGeometry(_parametric(triangle), 2)

From de2c784b95bea80e9851fdf44ec226d3ab6cff63 Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 22:17:41 -0500
Subject: [PATCH 29/33] Remove obsolete derivation

---
 docs/src/specializations.md | 40 -------------------------------------
 1 file changed, 40 deletions(-)

diff --git a/docs/src/specializations.md b/docs/src/specializations.md
index dc873601..ab0da442 100644
--- a/docs/src/specializations.md
+++ b/docs/src/specializations.md
@@ -11,43 +11,3 @@ There are several notable exceptions to how Meshes.jl defines [parametric functi
 - `Meshes.Rope` is a composite type that lacks a parametric function, but can be decomposed into `Meshes.Segment`s and then integrated by adding together the individual integrals.
 - `Meshes.Triangle` has a parametric function that takes coordinates on a 2D barycentric coordinate system. So, for `(::Meshes.Triangle)(t1, t2)`, the coordinates must obey: $t_1, t_2 \in [0,1]$ where $t_1 + t_2 \le 1$.
 - `Meshes.Tetrahedron` has a parametric function that takes coordinates on a 3D barycentric coordinate system. So, for `(::Meshes.Tetrahedron)(t1, t2)`, the coordinates must obey: $t_1, t_2, t_3 \in [0,1]$ where $t_1 + t_2 + t_3 \le 1$.
-
-## Triangle
-
-For a specified `Meshes.Triangle` surface with area $A$, let $u$ and $v$ be Barycentric coordinates that span the surface.
-```math
-\int_\triangle f(\bar{r}) \, \text{d}A
-    = \iint_\triangle f\left( \bar{r}(u,v) \right) \, \left( \text{d}u \wedge \text{d}v \right)
-```
-
-Since the geometric transformation from the originally-arbitrary domain to a Barycentric domain is linear, the magnitude of the surface element $\text{d}u \wedge \text{d}v$ is constant throughout the integration domain. This constant will be equal to twice the magnitude of $A$.
-```math
-\int_\triangle f(\bar{r}) \, \text{d}A
-    = 2A \int_0^1 \int_0^{1-v} f\left( \bar{r}(u,v) \right) \, \text{d}u \, \text{d}v
-```
-
-Since the integral domain is a right-triangle in the Barycentric domain, a nested application of Gauss-Kronrod quadrature rules is capable of computing the result, albeit inefficiently. However, many numerical integration methods that require rectangular bounds can not be directly applied.
-
-In order to enable integration methods that operate over rectangular bounds, two coordinate system transformations are applied: the first maps from Barycentric coordinates $(u, v)$ to polar coordinates $(r, \phi)$, and the second is a non-linear map from polar coordinates to a new curvilinear basis $(R, \phi)$.
-
-For the first transformation, let $u = r~\cos\phi$ and $v = r~\sin\phi$ where $\text{d}u~\text{d}v = r~\text{d}r~\text{d}\phi$. The Barycentric triangle's hypotenuse boundary line is described by the function $v(u) = 1 - u$. Substituting in the previous definitions leads to a new boundary line function in polar coordinate space $r(\phi) = 1 / (\sin\phi + \cos\phi)$.
-```math
-\int_0^1 \int_0^{1-v} f\left( \bar{r}(u,v) \right) \, \text{d}u \, \text{d}v =
-    \int_0^{\pi/2} \int_0^{1/(\sin\phi+\cos\phi)} f\left( \bar{r}(r,\phi) \right) \, r \, \text{d}r \, \text{d}\phi
-```
-
-These integral boundaries remain non-rectangular, so an additional transformation will be applied to a curvilinear $(R, \phi)$ space that normalizes all of the hypotenuse boundary line points to $R=1$. To achieve this, a function $R(r,\phi)$ is required such that $R(r_0, \phi) = 1$ where $r_0 = 1 / (\sin\phi + \cos\phi)$
-
-To achieve this, let $R(r, \phi) = r~(\sin\phi + \cos\phi)$. Now, substituting some terms leads to
-```math
-\int_0^{\pi/2} \int_0^{1/(\sin\phi+\cos\phi)} f\left( \bar{r}(r,\phi) \right) \, r \, \text{d}r \, \text{d}\phi
-    = \int_0^{\pi/2} \int_0^{r_0} f\left( \bar{r}(r,\phi) \right) \, \left(\frac{R}{\sin\phi + \cos\phi}\right) \, \text{d}r \, \text{d}\phi
-```
-
-Since $\text{d}R/\text{d}r = \sin\phi + \cos\phi$, a change of integral domain leads to
-```math
-\int_0^{\pi/2} \int_0^{r_0} f\left( \bar{r}(r,\phi) \right) \, \left(\frac{R}{\sin\phi + \cos\phi}\right) \, \text{d}r \, \text{d}\phi
-    = \int_0^{\pi/2} \int_0^1 f\left( \bar{r}(R,\phi) \right) \, \left(\frac{R}{\left(\sin\phi + \cos\phi\right)^2}\right) \, \text{d}R \, \text{d}\phi
-```
-
-The second term in this new integrand function serves as a correction factor that corrects for the impact of the non-linear domain transformation. Since all of the integration bounds are now constants, specialized integration methods can be defined for triangles that performs these domain transformations and then solve the new rectangular integration problem using a wider range of solver options.

From 4b3ad9dd021f49d0e977ea846b5edc03fa5e644d Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Sun, 1 Dec 2024 22:18:07 -0500
Subject: [PATCH 30/33] Update src/specializations/Tetrahedron.jl

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 src/specializations/Tetrahedron.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl
index 1aedfaf5..f074a2fd 100644
--- a/src/specializations/Tetrahedron.jl
+++ b/src/specializations/Tetrahedron.jl
@@ -30,7 +30,7 @@ end
 function _parametric(tetrahedron::Meshes.Tetrahedron)
     function f(t1, t2, t3)
         if any(Iterators.map(n -> (n < 0) || (n > 1), (t1, t2, t3)))
-            msg = "tetrahedron(t1, t2, t3) is not defined for t outside [0, 1]."
+            msg = "tetrahedron(t1, t2, t3) is not defined for (t1, t2, t3) outside [0, 1]³."
             throw(DomainError((t1, t2, t3), msg))
         end
 

From c3637cd2a82cb1de183b6934cf521729745ec48f Mon Sep 17 00:00:00 2001
From: Mike Ingold <mike.ingold@gmail.com>
Date: Mon, 2 Dec 2024 12:54:27 -0500
Subject: [PATCH 31/33] Remove docstring ref to Analytical

---
 src/differentiation.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/differentiation.jl b/src/differentiation.jl
index 5a5a0a6a..95232357 100644
--- a/src/differentiation.jl
+++ b/src/differentiation.jl
@@ -9,7 +9,7 @@ A category of types used to specify the desired method for calculating derivativ
 Derivatives are used to form Jacobian matrices when calculating the differential
 element size throughout the integration region.
 
-See also [`FiniteDifference`](@ref), [`Analytical`](@ref).
+See also [`FiniteDifference`](@ref).
 """
 abstract type DifferentiationMethod end
 

From 5b313f5229405e2124dc5ab187b09f929bb1447e Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Mon, 2 Dec 2024 13:23:58 -0500
Subject: [PATCH 32/33] Update src/specializations/Triangle.jl

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 src/specializations/Triangle.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index 04308ed8..22857df9 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -34,7 +34,8 @@ function _parametric(triangle::Meshes.Triangle)
             throw(DomainError((t1, t2), msg))
         end
 
-        return triangle(t1 * t2, t2 - (t1 * t2))
+        t1t2 = t1 * t2
+        return triangle(t1t2, t2 - t1t2)
     end
     return f
 end

From 29346fa0b2c42942164c00f9e264471f51cc9c17 Mon Sep 17 00:00:00 2001
From: Michael Ingold <mike.ingold@gmail.com>
Date: Mon, 2 Dec 2024 13:24:06 -0500
Subject: [PATCH 33/33] Update src/specializations/Triangle.jl

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 src/specializations/Triangle.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index 22857df9..26f3ea15 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -16,7 +16,7 @@ function integral(
         kwargs...
 )
     # Generate a _ParametricGeometry whose parametric function domain spans [0,1]²
-    param_triangle = _ParametricGeometry(_parametric(triangle), 2)
+    param_triangle = _ParametricGeometry(_parametric(triangle), Meshes.paramdim(triangle))
 
     # Integrate the _ParametricGeometry using the standard methods
     return _integral(f, param_triangle, rule; kwargs...)