Generalize differential element calculation to n-dims (#101)

* Add CliffordNumbers to deps

* Add _clifford utility

* Bugfix

* Change _clifford to _kvector, avoid an allocation

* Explicit namespace

* Add _Vec utility

* Generalize differential to n-dims

* using CliffordNumbers

* Explicit namespace, disable _Vec utility

* Bugfix

* Switch CliffordNumbers from import to using

* Remove _Vec utility

* Drop stale explicit import

* Convert LinearAlgebra to an import

* Enable full testing of 4D Box

* Analytic integrand/solution for 4D Box

* Bugfix

* Error when using GK rule >3D

* Bugfix

* Add explicit rtol for 4D Box

* Bump minimum CliffordNumbers (downgrade failing)

* Apply suggestions from code review

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

* Bump minimum again

* Convert all generic errors to specific errors

* Apply suggestions from code review

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

* Try new method for units handling

* Bugfix

* Generalize _units to all geometry types

* Condense generic GaussKronrod cases above 2D

* Apply format suggestion [skip ci]

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

* Update obsoleted comment

* Add tests for utils

* Change jacobian API to return via ntuple

* Apply format suggestions

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

* Remove _ukvector [skip ci]

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

* Organization

* Remove tests for _ukvector

* Remove @inline

* Increase minimum CliffordNumbers version

---------

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

Author: 3 people

Project.toml

@@ -3,6 +3,7 @@ uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
 version = "0.14.1"
 
 [deps]
+CliffordNumbers = "3998ac73-6bd4-4031-8035-f167dd3ed523"
 CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
@@ -12,6 +13,7 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
+CliffordNumbers = "0.1.4"
 CoordRefSystems = "0.12, 0.13, 0.14, 0.15"
 FastGaussQuadrature = "1"
 HCubature = "1.5"

src/MeshIntegrals.jl

@@ -1,9 +1,10 @@
 module MeshIntegrals
+using CliffordNumbers: CliffordNumbers, VGA, ∧
 using CoordRefSystems: CoordRefSystems, CRS
-using LinearAlgebra: LinearAlgebra, norm, ×, ⋅
 
 import FastGaussQuadrature
 import HCubature
+import LinearAlgebra
 import Meshes
 import QuadGK
 import Unitful

src/differentiation.jl

@@ -45,7 +45,7 @@ function jacobian(
         end
     end
 
-    return map(n -> ∂ₙr(ts, n), 1:Dim)
+    return ntuple(n -> ∂ₙr(ts, n), Dim)
 end
 
 function jacobian(
@@ -75,7 +75,7 @@ function jacobian(
     sigma(i) = B(i, N - 1)(t) .* (P[(i + 1) + 1] - P[(i) + 1])
     derivative = N .* sum(sigma, 0:(N - 1))
 
-    return [derivative]
+    return (derivative,)
 end
 
 ################################################################################
@@ -91,17 +91,16 @@ function for `geometry` at arguments `ts`.
 function differential(
         geometry::G,
         ts::V
-) where {G <: Meshes.Geometry, V <: Union{AbstractVector, Tuple}}
+) where {M, CRS, G <: Meshes.Geometry{M, CRS}, V <: Union{AbstractVector, Tuple}}
+    # Calculate the Jacobian, convert Vec -> KVector
     J = jacobian(geometry, ts)
+    J_kvecs = Iterators.map(_kvector, J)
 
-    # 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
+    # Extract units from Geometry type
+    Dim = Meshes.paramdim(geometry)
+    units = _units(geometry)^Dim
+
+    # Return norm of the exterior products
+    element = foldl(∧, J_kvecs)
+    return LinearAlgebra.norm(element) * units
 end

src/integral.jl

@@ -48,8 +48,9 @@ function _integral(
         return _integral_gk_1d(f, geometry, rule; kwargs...)
     elseif N == 2
         return _integral_gk_2d(f, geometry, rule; kwargs...)
-    elseif N == 3
-        return _integral_gk_3d(f, geometry, rule; kwargs...)
+    else
+        _error_unsupported_combination("geometry with more than two parametric dimensions",
+            "GaussKronrod")
     end
 end
 
@@ -126,13 +127,3 @@ function _integral_gk_2d(
     ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1]
     return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1]
 end
-
-# Integrating volumes with GaussKronrod not supported by default
-function _integral_gk_3d(
-        f,
-        geometry,
-        rule::GaussKronrod;
-        FP::Type{T} = Float64
-) where {T <: AbstractFloat}
-    error("Integrating this volume type with GaussKronrod not supported.")
-end

src/integral_aliases.jl

@@ -25,7 +25,8 @@ function lineintegral(
     if N == 1
         return integral(f, geometry, rule; kwargs...)
     else
-        error("Performing a line integral on a geometry with $N parametric dimensions not supported.")
+        throw(ArgumentError("Performing a line integral on a geometry \
+                            with $N parametric dimensions not supported."))
     end
 end
 
@@ -56,7 +57,8 @@ function surfaceintegral(
     if N == 2
         return integral(f, geometry, rule; kwargs...)
     else
-        error("Performing a surface integral on a geometry with $N parametric dimensions not supported.")
+        throw(ArgumentError("Performing a surface integral on a geometry \
+                            with $N parametric dimensions not supported."))
     end
 end
 
@@ -87,6 +89,7 @@ function volumeintegral(
     if N == 3
         return integral(f, geometry, rule; kwargs...)
     else
-        error("Performing a volume integral on a geometry with $N parametric dimensions not supported.")
+        throw(ArgumentError("Performing a volume integral on a geometry \
+                            with $N parametric dimensions not supported."))
     end
 end

src/specializations/Tetrahedron.jl

@@ -15,7 +15,7 @@ function integral(
         rule::GaussLegendre;
         FP::Type{T} = Float64
 ) where {F <: Function, T <: AbstractFloat}
-    error("Integrating a Tetrahedron with GaussLegendre not supported.")
+    _error_unsupported_combination("Tetrahedron", "GaussLegendre")
 end
 
 function integral(
@@ -40,5 +40,5 @@ function integral(
         rule::HAdaptiveCubature;
         FP::Type{T} = Float64
 ) where {F <: Function, T <: AbstractFloat}
-    error("Integrating a Tetrahedron with HAdaptiveCubature not supported.")
+    _error_unsupported_combination("Tetrahedron", "HAdaptiveCubature")
 end

src/utils.jl

@@ -8,5 +8,23 @@ function _gausslegendre(T, n)
     return T.(xs), T.(ws)
 end
 
-# Extract the length units used by the CRS of a Point
-_units(pt::Meshes.Point{M, CRS}) where {M, CRS} = first(CoordRefSystems.units(CRS))
+# Extract the length units used by the CRS of a Geometry
+function _units(g::Meshes.Geometry{M, CRS}) where {M, CRS}
+    return Unitful.unit(CoordRefSystems.lentype(CRS))
+end
+
+# Common error message structure
+function _error_unsupported_combination(geometry, rule)
+    msg = "Integrating a $geometry using a $rule rule not supported."
+    throw(ArgumentError(msg))
+end
+
+################################################################################
+#                        CliffordNumbers Interface
+################################################################################
+
+# Meshes.Vec -> ::CliffordNumber.KVector
+function _kvector(v::Meshes.Vec{Dim, T}) where {Dim, T}
+    ucoords = Iterators.map(Unitful.ustrip, v.coords)
+    return CliffordNumbers.KVector{1, VGA(Dim)}(ucoords...)
+end

test/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -12,6 +13,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 Aqua = "0.7, 0.8"
 ExplicitImports = "1.6.0"
+LinearAlgebra = "1"
 Meshes = "0.50, 0.51"
 QuadGK = "2.1.1"
 SpecialFunctions = "2"

test/combinations.jl

@@ -169,14 +169,32 @@ end
 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))
+    a = π
+    box = Box(Point(0, 0, 0, 0), Point(a, a, a, a))
 
-    f = p -> one(Float64)
+    function f(p::P) where {P <: Meshes.Point}
+        x1, x2, x3, x4 = ustrip.(to(p).coords)
+        σ(x) = sqrt(a^2 - x^2)
+        (σ(x1) + σ(x2) + σ(x3) + σ(x4)) * u"Ω/m^4"
+    end
+    fv(p) = fill(f(p), 3)
 
-    # Test for currently-unsupported >3D differentials
-    @test integral(f, box)≈1.0u"m^4" broken=true
+    # Scalar integrand
+    sol = 4a^3 * (π * 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(rtol = 1e-6))≈sol rtol=1e-6
+
+    # 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(rtol = 1e-6))≈vsol rtol=1e-6
+
+    # Integral aliases
+    @test_throws "not supported" lineintegral(f, box)
+    @test_throws "not supported" surfaceintegral(f, box)
+    @test_throws "not supported" volumeintegral(f, box)
 
     # Test jacobian with wrong number of parametric coordinates
     @test_throws ArgumentError jacobian(box, zeros(2))

test/utils.jl

@@ -0,0 +1,11 @@
+@testitem "Utilities" setup=[Setup] begin
+    using LinearAlgebra: norm
+
+    # _kvector
+    v = Meshes.Vec(3, 4)
+    @test norm(MeshIntegrals._kvector(v)) ≈ 5.0
+
+    # _units
+    p = Point(1.0u"cm", 2.0u"mm", 3.0u"m")
+    @test MeshIntegrals._units(p) == u"m"
+end