Rewrite of jacobian approximation (#96)

* Add tests for 1D Box

* Add units

* Adjust rtol

* Make auto tests verbose

* Add tests for 2D Box

* Apply format suggestions [skip ci]

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

* Tupleify

* Add tests for 3D Box

* Condense functions

* Remove reusable buffer

* Remove unsupported GK rule on 3D Box

* Convert ev to a generated ntuple

* Implement non-allocating jacobian

* Add explicit rtol

* Simplify Box constructors

* Remove old Box tests and obsolete lines

* Add dimension check and fix type stability issue

* Whitespace fixes [skip ci]

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

* Explicit namespace

* Remove broken status

* Add test for jacobian error condition

* Cleaner abstraction/iterators

* Apply suggestions from code review

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

* Fix 4D Box test

* Figured out destructuring of anonymous args

---------

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

Author: 3 people

src/differentiation.jl

@@ -18,49 +18,34 @@ function jacobian(
         ts::V;
         ε = 1e-6
 ) where {G <: Meshes.Geometry, V <: Union{AbstractVector, Tuple}}
-    T = eltype(ts)
+    Dim = Meshes.paramdim(geometry)
+    if Dim != length(ts)
+        throw(ArgumentError("ts must have same number of dimensions as geometry."))
+    end
 
-    # Get the partial derivative along the n'th axis via finite difference approximation
-    #   where ts is the current parametric position (εv is a reusable buffer)
-    function ∂ₙr!(εv, ts, n)
+    T = eltype(ts)
+    ε = T(ε)
+
+    # Get the partial derivative along the n'th axis via finite difference
+    #   approximation, where ts is the current parametric position
+    function ∂ₙr(ts, n)
+        # Build left/right parametric coordinates with non-allocating iterators 
+        left = Iterators.map(((i, t),) -> i == n ? t - ε : t, enumerate(ts))
+        right = Iterators.map(((i, t),) -> i == n ? t + ε : t, enumerate(ts))
+        # Select orientation of finite-diff
         if ts[n] < T(0.01)
-            return ∂ₙr_right!(εv, ts, n)
+            # Right
+            return (geometry(right...) - geometry(ts...)) / ε
         elseif T(0.99) < ts[n]
-            return ∂ₙr_left!(εv, ts, n)
+            # Left
+            return (geometry(ts...) - geometry(left...)) / ε
         else
-            return ∂ₙr_central!(εv, ts, n)
+            # Central
+            return (geometry(right...) - geometry(left...)) / 2ε
         end
     end
 
-    # Central finite difference
-    function ∂ₙr_central!(εv, ts, n)
-        εv[n] = ε
-        a = ts .- εv
-        b = ts .+ εv
-        return (geometry(b...) - geometry(a...)) / 2ε
-    end
-
-    # Left finite difference
-    function ∂ₙr_left!(εv, ts, n)
-        εv[n] = ε
-        a = ts .- εv
-        b = ts
-        return (geometry(b...) - geometry(a...)) / ε
-    end
-
-    # Right finite difference
-    function ∂ₙr_right!(εv, ts, n)
-        εv[n] = ε
-        a = ts
-        b = ts .+ εv
-        return (geometry(b...) - geometry(a...)) / ε
-    end
-
-    # Allocate a re-usable ε vector
-    εv = zeros(T, length(ts))
-
-    ∂ₙr(n) = ∂ₙr!(εv, ts, n)
-    return map(∂ₙr, 1:length(ts))
+    return map(n -> ∂ₙr(ts, n), 1:Dim)
 end
 
 function jacobian(

test/auto_tests.jl

@@ -39,8 +39,6 @@
 
     # Generate a @testset for item
     function autotest(item::SupportItem)
-        #@assert item.type == coordtype(item.geometry) "Item type mismatch"
-
         N = (item.type == Float32) ? 1000 : 100
         algorithm_set = [
             (GaussLegendre(N), item.gausslegendre),
@@ -77,15 +75,10 @@ end
     pt_n(T) = Point(T(0), T(1), T(0))
     pt_w(T) = Point(T(-1), T(0), T(0))
     pt_e(T) = Point(T(1), T(0), T(0))
-    pt_s(T) = Point(T(0), T(-1), T(0))
-    pt_z(T) = Point(T(0), T(0), T(1))
 
     # Test Geometries
     ball2d(T) = Ball(origin2d(T), T(2.0))
     ball3d(T) = Ball(origin3d(T), T(2.0))
-    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)))
     circle(T) = Circle(plane_xy(T), T(2.5))
     cyl(T) = Cylinder(pt_e(T), pt_w(T), T(2.5))
     cylsurf(T) = CylinderSurface(pt_e(T), pt_w(T), T(2.5))
@@ -101,24 +94,19 @@ 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("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),
         SupportItem("Circle{$T}", T, circle(T), 1, 1, 0, 0, 1, 1, 1),
         SupportItem("Cylinder{$T}", T, cyl(T), 1, 0, 0, 1, 1, 0, 1),
         SupportItem("CylinderSurface{$T}", T, cylsurf(T), 1, 0, 1, 0, 1, 1, 1),
         SupportItem("Disk{$T}", T, disk(T), 1, 0, 1, 0, 1, 1, 1),
-        # Frustum -- not yet supported
         SupportItem("ParaboloidSurface{$T}", T, parab(T), 1, 0, 1, 0, 1, 1, 1),
-        # SimpleMesh
         SupportItem("Sphere{2,$T}", T, sphere2d(T), 1, 1, 0, 0, 1, 1, 1),
         SupportItem("Sphere{3,$T}", T, sphere3d(T), 1, 0, 1, 0, 1, 1, 1),
         SupportItem("Tetrahedron", T, tetra(T), 1, 0, 0, 1, 0, 1, 0),
         SupportItem("Triangle{$T}", T, triangle(T), 1, 0, 1, 0, 1, 1, 1),
         SupportItem("Torus{$T}", T, torus(T), 1, 0, 1, 0, 1, 1, 1)
     ]
 
-    @testset "Float64 Geometries" begin
+    @testset "Float64 Geometries" verbose=true begin
         map(autotest, SUPPORT_MATRIX(Float64))
     end
 end

test/combinations.jl

@@ -34,10 +34,102 @@
     @test_throws DomainError jacobian(curve, [1.1])
 end
 
-@testitem "Meshes.Box" setup=[Setup] begin
+@testitem "Meshes.Box 1D" setup=[Setup] begin
+    a = π
+    box = Box(Point(0), Point(a))
+
+    function f(p::P) where {P <: Meshes.Point}
+        t = ustrip(p.coords.x)
+        sqrt(a^2 - t^2) * u"Ω/m"
+    end
+    fv(p) = fill(f(p), 3)
+
+    # Scalar integrand
+    sol = π * a^2 / 4 * u"Ω"
+    @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6
+    @test integral(f, box, GaussKronrod()) ≈ sol
+    @test integral(f, box, HAdaptiveCubature()) ≈ sol
+
+    # Vector integrand
+    vsol = fill(sol, 3)
+    @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6
+    @test integral(fv, box, GaussKronrod()) ≈ vsol
+    @test integral(fv, box, HAdaptiveCubature()) ≈ vsol
+
+    # Integral aliases
+    @test lineintegral(f, box) ≈ sol
+    @test_throws "not supported" surfaceintegral(f, box)
+    @test_throws "not supported" volumeintegral(f, box)
+end
+
+@testitem "Meshes.Box 2D" setup=[Setup] begin
+    a = π
+    box = Box(Point(0, 0), Point(a, a))
+
+    function f(p::P) where {P <: Meshes.Point}
+        x, y = ustrip.((p.coords.x, p.coords.y))
+        (sqrt(a^2 - x^2) + sqrt(a^2 - y^2)) * u"Ω/m^2"
+    end
+    fv(p) = fill(f(p), 3)
+
+    # Scalar integrand
+    sol = 2a * (π * a^2 / 4) * u"Ω"
+    @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6
+    @test integral(f, box, GaussKronrod()) ≈ sol
+    @test integral(f, box, HAdaptiveCubature()) ≈ sol
+
+    # Vector integrand
+    vsol = fill(sol, 3)
+    @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6
+    @test integral(fv, box, GaussKronrod()) ≈ vsol
+    @test integral(fv, box, HAdaptiveCubature()) ≈ vsol
+
+    # Integral aliases
+    @test_throws "not supported" lineintegral(f, box)
+    @test surfaceintegral(f, box) ≈ sol
+    @test_throws "not supported" volumeintegral(f, box)
+end
+
+@testitem "Meshes.Box 3D" setup=[Setup] begin
+    a = π
+    box = Box(Point(0, 0, 0), Point(a, a, a))
+
+    function f(p::P) where {P <: Meshes.Point}
+        x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z))
+        (sqrt(a^2 - x^2) + sqrt(a^2 - y^2) + sqrt(a^2 - z^2)) * u"Ω/m^3"
+    end
+    fv(p) = fill(f(p), 3)
+
+    # Scalar integrand
+    sol = 3a^2 * (π * a^2 / 4) * u"Ω"
+    @test integral(f, box, GaussLegendre(100))≈sol rtol=1e-6
+    @test_throws "not supported" integral(f, box, GaussKronrod())
+    @test integral(f, box, HAdaptiveCubature()) ≈ sol
+
+    # Vector integrand
+    vsol = fill(sol, 3)
+    @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6
+    @test_throws "not supported" integral(fv, box, GaussKronrod())
+    @test integral(fv, box, HAdaptiveCubature()) ≈ vsol
+
+    # Integral aliases
+    @test_throws "not supported" lineintegral(f, box)
+    @test_throws "not supported" surfaceintegral(f, box)
+    @test volumeintegral(f, box) ≈ sol
+end
+
+@testitem "Meshes.Box 4D" setup=[Setup] begin
+    a = zero(Float64)
+    b = zero(Float64)
+    box = Box(Point(a, a, a, a), Point(b, b, b, b))
+
+    f = p -> one(Float64)
+
     # Test for currently-unsupported >3D differentials
-    box4d = Box(Point(zeros(4)...), Point(ones(4)...))
-    @test integral(f -> one(Float64), box4d)≈1.0u"m^4" broken=true
+    @test integral(f, box)≈1.0u"m^4" broken=true
+
+    # Test jacobian with wrong number of parametric coordinates
+    @test_throws ArgumentError jacobian(box, zeros(2))
 end
 
 @testitem "Meshes.Cone" setup=[Setup] begin

test/floating_point_types.jl

@@ -11,33 +11,33 @@ end
 
 @testitem "Alternate floating types" setup=[Setup, BaseAtol] begin
     @testset "$FP" for FP in (Float32, BigFloat)
-        # typeof @test's are currently broken for Float32, see GitHub Issue 74
-
         # Rectangular volume with unit integrand
         f = p -> one(FP)
-        box1d = Box(Point(fill(zero(FP) * u"m", 1)...), Point(fill(one(FP) * u"m", 1)...))
-        box2d = Box(Point(fill(zero(FP) * u"m", 2)...), Point(fill(one(FP) * u"m", 2)...))
-        box3d = Box(Point(fill(zero(FP) * u"m", 3)...), Point(fill(one(FP) * u"m", 3)...))
+        a = zero(FP) * u"m"
+        b = one(FP) * u"m"
+        box1d = Box(Point(a), Point(b))
+        box2d = Box(Point(a, a), Point(b, b))
+        box3d = Box(Point(a, a, a), Point(b, b, b))
 
         # Check HCubature integrals (same method invoked for all dimensions)
         int_HC = integral(f, box1d, HAdaptiveCubature(); FP = FP)
         @test int_HC≈one(FP) * u"m" atol=baseatol[FP] * u"m"
-        @test typeof(int_HC.val)==FP broken=(FP == Float32)
+        @test typeof(int_HC.val) == FP
 
         # Check Gauss-Legendre integral in 1D
         int_GL_1D = integral(f, box1d, GaussLegendre(100); FP = FP)
         @test int_GL_1D≈one(FP) * u"m" atol=baseatol[FP] * u"m"
-        @test typeof(int_GL_1D.val)==FP broken=(FP == Float32)
+        @test typeof(int_GL_1D.val) == FP
 
         # Check Gauss-Legendre integral in 2D
         int_GL_2D = integral(f, box2d, GaussLegendre(100); FP = FP)
         @test int_GL_2D≈one(FP) * u"m^2" atol=2baseatol[FP] * u"m^2"
-        @test typeof(int_GL_2D.val)==FP broken=(FP == Float32)
+        @test typeof(int_GL_2D.val) == FP
 
         # Check Gauss-Legendre integral in 3D
         int_GL_3D = integral(f, box3d, GaussLegendre(100); FP = FP)
         @test int_GL_3D≈one(FP) * u"m^3" atol=3baseatol[FP] * u"m^3"
-        @test typeof(int_GL_3D.val)==FP broken=(FP == Float32)
+        @test typeof(int_GL_3D.val) == FP
     end
 end