Implement differential and jacobian for Bezier curve integrals (#92)

* Convert Bezier methods from static dl to using derivative

* apply format suggestions

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

* fix dimension of test point

* Add new tests for BezierCurve

* Merge derivative with jacobian

* Code formatting

* Update names

* Bugfix - deconstruct t from ts [skip ci]

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

* Apply format suggestions [skip ci]

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

* Clarify variable name

* Add argument type parameters to differential

* Bugfix - missing comma to capture as tuple

* Updated path name

* Remove old tests

* Add an rtol for tests

* Update rtol values

* Update rtol values

* Add jacobian to API page

---------

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

Author: 3 people

docs/src/api.md

@@ -21,3 +21,9 @@ MeshIntegrals.GaussKronrod
 MeshIntegrals.GaussLegendre
 MeshIntegrals.HAdaptiveCubature
 ```
+
+## Derivatives
+
+```@docs
+MeshIntegrals.jacobian
+```

src/differentiation.jl

@@ -14,10 +14,10 @@ finite-difference approximation with step size `ε`.
 - `ε`: step size to use for the finite-difference approximation
 """
 function jacobian(
-        geometry,
-        ts;
+        geometry::G,
+        ts::V;
         ε = 1e-6
-)
+) where {G <: Meshes.Geometry, V <: Union{AbstractVector, Tuple}}
     T = eltype(ts)
 
     # Get the partial derivative along the n'th axis via finite difference approximation
@@ -63,44 +63,11 @@ function jacobian(
     return map(∂ₙr, 1:length(ts))
 end
 
-"""
-    differential(geometry, ts)
-
-Calculate the differential element (length, area, volume, etc) of the parametric
-function for `geometry` at arguments `ts`.
-"""
-function differential(
-        geometry,
-        ts
-)
-    J = jacobian(geometry, ts)
-
-    # TODO generalize this with geometric algebra, e.g.: norm(foldl(∧, J))
-    if length(J) == 1
-        return norm(J[1])
-    elseif length(J) == 2
-        return norm(J[1] × J[2])
-    elseif length(J) == 3
-        return abs((J[1] × J[2]) ⋅ J[3])
-    else
-        error("Not implemented yet. Please report this as an Issue on GitHub.")
-    end
-end
-
-################################################################################
-#                             Partial Derivatives
-################################################################################
-
-"""
-    derivative(b::BezierCurve, t)
-
-Determine the vector derivative of a Bezier curve `b` for the point on the
-curve parameterized by value `t`.
-"""
-function derivative(
+function jacobian(
         bz::Meshes.BezierCurve,
-        t
-)
+        ts::V
+) where {V <: Union{AbstractVector, Tuple}}
+    t = only(ts)
     # Parameter t restricted to domain [0,1] by definition
     if t < 0 || t > 1
         throw(DomainError(t, "b(t) is not defined for t outside [0, 1]."))
@@ -121,5 +88,35 @@ function derivative(
     # Derivative = N Σ_{i=0}^{N-1} sigma(i)
     #   P indices adjusted for Julia 1-based array indexing
     sigma(i) = B(i, N - 1)(t) .* (P[(i + 1) + 1] - P[(i) + 1])
-    return N .* sum(sigma, 0:(N - 1))
+    derivative = N .* sum(sigma, 0:(N - 1))
+
+    return [derivative]
+end
+
+################################################################################
+#                          Differential Elements
+################################################################################
+
+"""
+    differential(geometry, ts)
+
+Calculate the differential element (length, area, volume, etc) of the parametric
+function for `geometry` at arguments `ts`.
+"""
+function differential(
+        geometry::G,
+        ts::V
+) where {G <: Meshes.Geometry, V <: Union{AbstractVector, Tuple}}
+    J = jacobian(geometry, ts)
+
+    # TODO generalize this with geometric algebra, e.g.: norm(foldl(∧, J))
+    if length(J) == 1
+        return norm(J[1])
+    elseif length(J) == 2
+        return norm(J[1] × J[2])
+    elseif length(J) == 3
+        return abs((J[1] × J[2]) ⋅ J[3])
+    else
+        error("Not implemented yet. Please report this as an Issue on GitHub.")
+    end
 end

src/integral.jl

@@ -87,14 +87,14 @@ function _integral(
 ) where {T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
 
-    integrand(t) = f(geometry(t...)) * differential(geometry, t)
+    integrand(ts) = f(geometry(ts...)) * differential(geometry, ts)
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
     testpoint_parametriccoord = zeros(FP, N)
     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))
+    uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts))
     # Integrate only the unitless values
     value = HCubature.hcubature(uintegrand, zeros(FP, N), ones(FP, N); rule.kwargs...)[1]
 
@@ -112,7 +112,7 @@ function _integral_gk_1d(
         rule::GaussKronrod;
         FP::Type{T} = Float64
 ) where {T <: AbstractFloat}
-    integrand(t) = f(geometry(t)) * differential(geometry, (t))
+    integrand(t) = f(geometry(t)) * differential(geometry, (t,))
     return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
 end
 

src/specializations/BezierCurve.jl

@@ -33,9 +33,10 @@ function integral(
     # Change of variables: x [-1,1] ↦ t [0,1]
     t(x) = FP(1 // 2) * x + FP(1 // 2)
     point(x) = curve(t(x), alg)
+    integrand(x) = f(point(x)) * differential(curve, (t(x),))
 
     # Integrate f along curve and apply domain-correction for [-1,1] ↦ [0, length]
-    return FP(1 // 2) * length(curve) * sum(w .* f(point(x)) for (w, x) in zip(ws, xs))
+    return FP(1 // 2) * sum(w .* integrand(x) for (w, x) in zip(ws, xs))
 end
 
 """
@@ -56,9 +57,8 @@ function integral(
         FP::Type{T} = Float64,
         alg::Meshes.BezierEvalMethod = Meshes.Horner()
 ) where {F <: Function, T <: AbstractFloat}
-    len = length(curve)
     point(t) = curve(t, alg)
-    integrand(t) = len * f(point(t))
+    integrand(t) = f(point(t)) * differential(curve, (t,))
     return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
 end
 
@@ -80,16 +80,15 @@ function integral(
         FP::Type{T} = Float64,
         alg::Meshes.BezierEvalMethod = Meshes.Horner()
 ) where {F <: Function, T <: AbstractFloat}
-    len = length(curve)
     point(t) = curve(t, alg)
-    integrand(t) = len * f(point(t[1]))
+    integrand(ts) = f(point(only(ts))) * differential(curve, ts)
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
-    testpoint_parametriccoord = zeros(FP, 3)
+    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))
+    uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts))
     # Integrate only the unitless values
     value = HCubature.hcubature(uintegrand, zeros(FP, 1), ones(FP, 1); rule.kwargs...)[1]
 

test/auto_tests.jl

@@ -83,7 +83,6 @@ end
     # Test Geometries
     ball2d(T) = Ball(origin2d(T), T(2.0))
     ball3d(T) = Ball(origin3d(T), T(2.0))
-    bezier(T) = BezierCurve([Point(cos(t), sin(t), 0) for t in T(0):T(0.1):T(2π)])
     box1d(T) = Box(Point(T(-1)), Point(T(1)))
     box2d(T) = Box(Point(T(-1), T(-1)), Point(T(1), T(1)))
     box3d(T) = Box(Point(T(-1), T(-1), T(-1)), Point(T(1), T(1), T(-1)))
@@ -102,7 +101,6 @@ end
         # Name, T type, example,    integral,line,surface,volume,    GaussLegendre,GaussKronrod,HAdaptiveCubature
         SupportItem("Ball{2,$T}", T, ball2d(T), 1, 0, 1, 0, 1, 1, 1),
         SupportItem("Ball{3,$T}", T, ball3d(T), 1, 0, 0, 1, 1, 0, 1),
-        SupportItem("BezierCurve{$T}", T, bezier(T), 1, 1, 0, 0, 1, 1, 1),
         SupportItem("Box{1,$T}", T, box1d(T), 1, 1, 0, 0, 1, 1, 1),
         SupportItem("Box{2,$T}", T, box2d(T), 1, 0, 1, 0, 1, 1, 1),
         SupportItem("Box{3,$T}", T, box3d(T), 1, 0, 0, 1, 1, 0, 1),

test/combinations.jl

@@ -2,6 +2,33 @@
 # - All supported combinations of integral(f, ::Geometry, ::IntegrationAlgorithm) produce accurate results
 # - Invalid applications of integral aliases (e.g. lineintegral) produce a descriptive error
 
+@testitem "Meshes.BezierCurve" setup=[Setup] begin
+    curve = BezierCurve(
+        [Point(t * u"m", sin(t) * u"m", 0.0u"m") for t in range(-pi, pi, length = 361)]
+    )
+
+    f(x, y, z) = (1 / sqrt(1 + cos(x)^2)) * u"Ω/m"
+    f(p::Meshes.Point) = f(ustrip(to(p))...)
+    fv(p) = fill(f(p), 3)
+
+    # Scalar integrand
+    sol = 2π * u"Ω"
+    @test integral(f, curve, GaussLegendre(100))≈sol rtol=0.5e-2
+    @test integral(f, curve, GaussKronrod())≈sol rtol=0.5e-2
+    @test integral(f, curve, HAdaptiveCubature())≈sol rtol=0.5e-2
+
+    # Vector integrand
+    vsol = fill(sol, 3)
+    @test integral(fv, curve, GaussLegendre(100))≈vsol rtol=0.5e-2
+    @test integral(fv, curve, GaussKronrod())≈vsol rtol=0.5e-2
+    @test integral(fv, curve, HAdaptiveCubature())≈vsol rtol=0.5e-2
+
+    # Integral aliases
+    @test lineintegral(f, curve)≈sol rtol=0.5e-2
+    @test_throws "not supported" surfaceintegral(f, curve)
+    @test_throws "not supported" volumeintegral(f, curve)
+end
+
 @testitem "Meshes.Cone" setup=[Setup] begin
     r = 2.5u"m"
     h = 2.5u"m"