Implement more _ParametricGeometry transforms, remove Analytical (#139)

* Implement _ParametricGeometry for Line

* Implement _ParametricGeometry for Plane

* Implement _ParametricGeometry for Ray

* Update and reorganize analytical tests

* Add comments

* Use paramdim for clarity

* Test _parametric -> anonymous function

* Update _parametric API to return an anonymous function

* Remove obsolete N, enable Tetrahedron test

* Remove hanging comma

* Update field name

* Add analytical derivative in comment

* Update _parametric for Triangle and Tetrahedron

* Update tests

* put parametrized to .typos.toml

* Implement _ParametricGeometry for Triangle

* Temp disable some Analytical tests

* Temp disable more tests

* Test disabling bound check

* Restore bounds check, try new formulation

* Remove old commented-out code

* Add comments

* Remove obsolete tests

* Remove obsolete utils and update default diff_method

* Remove obsolete imports

* Remove obsolete Analytical

* Remove obsolete export

* Formatting

* Remove obsolete derivation

* Update src/specializations/Tetrahedron.jl

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

* Remove docstring ref to Analytical

* Update src/specializations/Triangle.jl

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

* Update src/specializations/Triangle.jl

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

---------

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

Author: mikeingold

.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
 

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.

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")
 

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
 
@@ -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

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